diff --git a/src/auxiliary/math.jl b/src/auxiliary/math.jl
index 9e3aaa181bf..23a4278add1 100644
--- a/src/auxiliary/math.jl
+++ b/src/auxiliary/math.jl
@@ -160,6 +160,19 @@ Given ε = 1.0e-4, we use the following algorithm.
     end
 end
 
+# TODO: This may need some performance tuning, e.g., to optimize the
+#       cutoff value `epsilon_f2` or the number of terms in the Taylor expansion
+@inline function ln_mean(x::Float32, y::Float32)
+    epsilon_f2 = 1.0f-4 # cut-off value to be changed for Float32
+    f2 = (x * (x - 2 * y) + y * y) / (x * (x + 2 * y) + y * y) # f2 = f^2
+    if f2 < epsilon_f2
+        # Taylor expansion to be changed for Float32
+        return (x + y) / @evalpoly(f2, 2.0f0, 2.0f0/3, 2.0f0/5, 2.0f0/7)
+    else
+        return (y - x) / log(y / x)
+    end
+end
+
 """
     inv_ln_mean(x, y)
 
@@ -180,6 +193,19 @@ multiplication.
     end
 end
 
+# TODO: This may need some performance tuning, e.g., to optimize the
+#       cutoff value `epsilon_f2` or the number of terms in the Taylor expansion
+@inline function inv_ln_mean(x::Float32, y::Float32)
+    epsilon_f2 = 1.0f-4 # cut-off value to be changed for Float32
+    f2 = (x * (x - 2 * y) + y * y) / (x * (x + 2 * y) + y * y) # f2 = f^2
+    if f2 < epsilon_f2
+        # Taylor expansion to be changed for Float32
+        return @evalpoly(f2, 2.0f0, 2.0f0/3, 2.0f0/5, 2.0f0/7) / (x + y)
+    else
+        return log(y / x) / (y - x)
+    end
+end
+
 # `Base.max` and `Base.min` perform additional checks for signed zeros and `NaN`s
 # which are not present in comparable functions in Fortran/C++. For example,
 # ```julia