Support and test fma, fma and fmna for complex number

Fix #387

Author: serge-sans-paille

include/xsimd/math/xsimd_math_complex.hpp

@@ -815,6 +815,34 @@ namespace xsimd
                 return hypot(z.real(), z.imag());
             }
 
+            static batch_type fma(const batch_type& a, const batch_type& b, const batch_type& c)
+            {
+                real_batch res_r = ::xsimd::fms(a.real(), b.real(), ::xsimd::fms(a.imag(), b.imag(), c.real()));
+                real_batch res_i = ::xsimd::fma(a.real(), b.imag(), ::xsimd::fma(a.imag(), b.real(), c.imag()));
+                return {res_r, res_i};
+            }
+
+            static batch_type fms(const batch_type& a, const batch_type& b, const batch_type& c)
+            {
+                real_batch res_r = ::xsimd::fms(a.real(), b.real(), ::xsimd::fma(a.imag(), b.imag(), c.real()));
+                real_batch res_i = ::xsimd::fma(a.real(), b.imag(), ::xsimd::fms(a.imag(), b.real(), c.imag()));
+                return {res_r, res_i};
+            }
+
+            static batch_type fnma(const batch_type& a, const batch_type& b, const batch_type& c)
+            {
+                real_batch res_r = - ::xsimd::fms(a.real(), b.real(), ::xsimd::fma(a.imag(), b.imag(), c.real()));
+                real_batch res_i = - ::xsimd::fma(a.real(), b.imag(), ::xsimd::fms(a.imag(), b.real(), c.imag()));
+                return {res_r, res_i};
+            }
+
+            static batch_type fnms(const batch_type& a, const batch_type& b, const batch_type& c)
+            {
+                real_batch res_r = - ::xsimd::fms(a.real(), b.real(), ::xsimd::fms(a.imag(), b.imag(), c.real()));
+                real_batch res_i = - ::xsimd::fma(a.real(), b.imag(), ::xsimd::fma(a.imag(), b.real(), c.imag()));
+                return {res_r, res_i};
+            }
+
             static batch_type sqrt(const batch_type& z)
             {
                 using rvt = typename real_batch::value_type;

include/xsimd/math/xsimd_scalar.hpp

@@ -375,6 +375,31 @@ namespace xsimd
         return std::fma(a, b, c);
     }
 
+    namespace detail
+    {
+        template <class C>
+        inline C fma_complex_scalar_impl(const C& a, const C& b, const C& c)
+        {
+            return {fms(a.real(), b.real(), fms(a.imag(), b.imag(), c.real())),
+                    fma(a.real(), b.imag(), fma(a.imag(), b.real(), c.imag()))};
+        }
+    }
+
+    template <class T>
+    inline std::complex<T> fma(const std::complex<T>& a, const std::complex<T>& b, const std::complex<T>& c)
+    {
+        return detail::fma_complex_scalar_impl(a, b, c);
+    }
+
+
+#ifdef XSIMD_ENABLE_XTL_COMPLEX
+    template <class T, bool i3ec>
+    inline xtl::xcomplex<T, T, i3ec> fma(const xtl::xcomplex<T, T, i3ec>& a, const xtl::xcomplex<T, T, i3ec>& b, const xtl::xcomplex<T, T, i3ec>& c)
+    {
+        return detail::fma_complex_scalar_impl(a, b, c);
+    }
+#endif
+
     inline void sincos(float val, float&s, float& c)
     {
 #if defined(__APPLE__)

test/test_batch_complex.cpp

@@ -468,6 +468,42 @@ class batch_complex_test : public testing::Test
         }
     }
 
+    void test_fused_operations() const
+    {
+        // fma
+        {
+            array_type expected;
+            std::transform(lhs.cbegin(), lhs.cend(), rhs.begin(), expected.begin(),
+                            [](const value_type& l, const value_type& r) { return l * r + r; });
+            batch_type res = xsimd::fma(batch_lhs(), batch_rhs(), batch_rhs());
+            EXPECT_BATCH_EQ(res, expected) << print_function_name("fma");
+        }
+        // fms
+        {
+            array_type expected;
+            std::transform(lhs.cbegin(), lhs.cend(), rhs.begin(), expected.begin(),
+                            [](const value_type& l, const value_type& r) { return l * r - r; });
+            batch_type res = fms(batch_lhs(), batch_rhs(), batch_rhs());
+            EXPECT_BATCH_EQ(res, expected) << print_function_name("fms");
+        }
+        // fnma
+        {
+            array_type expected;
+            std::transform(lhs.cbegin(), lhs.cend(), rhs.begin(), expected.begin(),
+                            [](const value_type& l, const value_type& r) { return -l * r + r; });
+            batch_type res = fnma(batch_lhs(), batch_rhs(), batch_rhs());
+            EXPECT_BATCH_EQ(res, expected) << print_function_name("fnma");
+        }
+        // fnms
+        {
+            array_type expected;
+            std::transform(lhs.cbegin(), lhs.cend(), rhs.begin(), expected.begin(),
+                            [](const value_type& l, const value_type& r) { return -l * r - r; });
+            batch_type res = fnms(batch_lhs(), batch_rhs(), batch_rhs());
+            EXPECT_BATCH_EQ(res, expected) << print_function_name("fnms");
+        }
+    }
+
 private:
 
     batch_type batch_lhs() const
@@ -522,3 +558,7 @@ TYPED_TEST(batch_complex_test, horizontal_operations)
     this->test_horizontal_operations();
 }
 
+TYPED_TEST(batch_complex_test, fused_operations)
+{
+    this->test_fused_operations();
+}