diff --git a/include/interpolation.h b/include/interpolation.h
index 5dd98ea3f..cbe274d42 100644
--- a/include/interpolation.h
+++ b/include/interpolation.h
@@ -32,6 +32,7 @@
 
 #include <math.h>
 #include "lmms_constants.h"
+#include "lmms_math.h"
 
 inline float hermiteInterpolate( float x0, float x1, float x2, float x3,
 								float frac_pos )
@@ -80,32 +81,13 @@ inline float cubicInterpolate( float v0, float v1, float v2, float v3, float x )
 inline float cosinusInterpolate( float v0, float v1, float x )
 {
 	const float f = ( 1.0f - cosf( x * F_PI ) ) * 0.5f;
-#ifdef FP_FAST_FMAF
-	#ifndef __clang__
-		return fmaf( f, v1-v0, v0 );
-	#else
-		return fma( f, v1-v0, v0 );
-	#endif
-#else
-	return f * (v1-v0) + v0;
-#endif	
-//	return( v0*f + v1*( 1.0f-f ) );
+	return fastFmaf( f, v1-v0, v0 );
 }
 
 
 inline float linearInterpolate( float v0, float v1, float x )
 {
-// take advantage of fma function if present in hardware
-
-#ifdef FP_FAST_FMAF
-	#ifndef __clang__
-		return fmaf( x, v1-v0, v0 );
-	#else
-		return fma( x, v1-v0, v0 );
-	#endif
-#else
-	return x * (v1-v0) + v0;
-#endif	
+	return fastFmaf( x, v1-v0, v0 );
 }
 
 
diff --git a/include/lmms_math.h b/include/lmms_math.h
index 47136ae94..3d23abcd3 100644
--- a/include/lmms_math.h
+++ b/include/lmms_math.h
@@ -130,7 +130,40 @@ static inline int fast_rand()
 	return( (unsigned)( next / 65536 ) % 32768 );
 }
 
+//! @brief Takes advantage of fmal() function if present in hardware
+static inline long double fastFmal( long double a, long double b, long double c ) {
+#ifdef FP_FAST_FMAL
+	#ifdef __clang__
+		return fma( a, b, c );
+	#else
+		return fmal( a, b, c );
+	#endif
+#else
+	return a * b + c;
+#endif
+}
 
+//! @brief Takes advantage of fmaf() function if present in hardware
+static inline float fastFmaf( float a, float b, float c ) {
+#ifdef FP_FAST_FMAF
+	#ifdef __clang__
+		return fma( a, b, c );
+	#else
+		return fmaf( a, b, c );
+	#endif
+#else
+	return a * b + c;
+#endif
+}
+
+//! @brief Takes advantage of fma() function if present in hardware
+static inline double fastFma( double a, double b, double c ) {
+#ifdef FP_FAST_FMA
+	return fma( a, b, c );
+#else
+	return a * b + c;
+#endif
+}
 
 // source: http://martin.ankerl.com/2007/10/04/optimized-pow-approximation-for-java-and-c-c/
 static inline double fastPow( double a, double b )