diff --git a/src/math.hpp b/src/math.hpp
index 4d8ee51c103cc74107938135b8094276fde10c47..bf964ba38c07487914ffcf5768b7f39aee63a5d1 100644
--- a/src/math.hpp
+++ b/src/math.hpp
@@ -1,9 +1,9 @@
 #pragma once
 
 #include <cmath>
+#include <limits>
 #include <utility>
 
-
 namespace nest {
 namespace mc {
 namespace math {
@@ -14,6 +14,12 @@ T constexpr pi()
     return T(3.1415926535897932384626433832795);
 }
 
+template <typename T = float>
+T constexpr infinity()
+{
+    return std::numeric_limits<T>::infinity();
+}
+
 template <typename T>
 T constexpr mean(T a, T b)
 {
diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt
index b6a7d811f24a88fdd618ea80e2e2ac17d6e919f5..60ca6ab0572c74a580399364609d16b9c242adaa 100644
--- a/tests/unit/CMakeLists.txt
+++ b/tests/unit/CMakeLists.txt
@@ -12,6 +12,7 @@ set(TEST_SOURCES
     test_cell_group.cpp
     test_lexcmp.cpp
     test_mask_stream.cpp
+    test_math.cpp
     test_matrix.cpp
     test_mechanisms.cpp
     test_nop.cpp
diff --git a/tests/unit/test_math.cpp b/tests/unit/test_math.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b32bd194cb75b9df1c31c34c0323fdc0b95bce9e
--- /dev/null
+++ b/tests/unit/test_math.cpp
@@ -0,0 +1,36 @@
+#include <cmath>
+#include <limits>
+
+#include "gtest.h"
+#include <math.hpp>
+
+using namespace nest::mc::math;
+
+TEST(math, infinity) {
+    // check values for float, double, long double
+    auto finf = infinity<float>();
+    EXPECT_TRUE((std::is_same<float, decltype(finf)>::value));
+    EXPECT_TRUE(std::isinf(finf));
+    EXPECT_GT(finf, 0.f);
+
+    auto dinf = infinity<double>();
+    EXPECT_TRUE((std::is_same<double, decltype(dinf)>::value));
+    EXPECT_TRUE(std::isinf(dinf));
+    EXPECT_GT(dinf, 0.0);
+
+    auto ldinf = infinity<long double>();
+    EXPECT_TRUE((std::is_same<long double, decltype(ldinf)>::value));
+    EXPECT_TRUE(std::isinf(ldinf));
+    EXPECT_GT(ldinf, 0.0l);
+
+    // check default value promotes correctly (i.e., acts like INFINITY)
+    struct {
+        float f;
+        double d;
+        long double ld;
+    } check = {infinity<>(), infinity<>(), infinity<>()};
+
+    EXPECT_EQ(std::numeric_limits<float>::infinity(), check.f);
+    EXPECT_EQ(std::numeric_limits<double>::infinity(), check.d);
+    EXPECT_EQ(std::numeric_limits<long double>::infinity(), check.ld);
+}