Support of other real types (Float32, auxiliary)

src/auxiliary/math.jl

@@ -160,6 +160,18 @@ Given ε = 1.0e-4, we use the following algorithm.
     end
 end
 
+# TODO: Performance tuning for `ln_mean` (Float32), please find below
+@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 +192,18 @@ multiplication.
     end
 end
 
+# TODO: Performance tuning for `inv_ln_mean` (Float32), please find below
+@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