C99: Initial fma support.

include/boost/multiprecision/detail/default_ops.hpp

@@ -385,7 +385,7 @@ inline void eval_multiply_default(T& t, const U& u, const V& v)
    }
    else
    {
-      t = u;
+      t = number<T>::canonical_value(u);
       eval_multiply(t, v);
    }
 }
@@ -395,13 +395,13 @@ inline void eval_multiply(T& t, const U& u, const V& v)
    eval_multiply_default(t, u, v);
 }
 
-template <class T, class U, class V, class X>
-inline typename disable_if_c<!is_same<T, U>::value && is_same<T, V>::value>::type eval_multiply_add(T& t, const U& u, const V& v, const X& x)
+template <class T>
+inline void eval_multiply_add(T& t, const T& u, const T& v, const T& x)
 {
    if((void*)&x == (void*)&t)
    {
       T z;
-      z = x;
+      z = number<T>::canonical_value(x);
       eval_multiply_add(t, u, v, z);
    }
    else
@@ -410,6 +410,25 @@ inline typename disable_if_c<!is_same<T, U>::value && is_same<T, V>::value>::typ
       eval_add(t, x);
    }
 }
+
+template <class T, class U>
+inline typename boost::disable_if_c<boost::is_same<T, U>::value, T>::type make_T(const U& u)
+{
+   T t;
+   t = number<T>::canonical_value(u);
+   return BOOST_MP_MOVE(t);
+}
+template <class T>
+inline const T& make_T(const T& t)
+{
+   return t;
+}
+
+template <class T, class U, class V, class X>
+inline typename disable_if_c<!is_same<T, U>::value && is_same<T, V>::value>::type eval_multiply_add(T& t, const U& u, const V& v, const X& x)
+{
+   eval_multiply_add(t, make_T<T>(u), make_T<T>(v), make_T<T>(x));
+}
 template <class T, class U, class V, class X>
 inline typename enable_if_c<!is_same<T, U>::value && is_same<T, V>::value>::type eval_multiply_add(T& t, const U& u, const V& v, const X& x)
 {
@@ -1000,7 +1019,7 @@ inline typename boost::enable_if_c<boost::is_arithmetic<A>::value>::type eval_fd
    typedef typename boost::multiprecision::detail::canonical<A, T>::type arithmetic_type;
    static const ui_type zero = 0u;
    arithmetic_type canonical_b = b;
-   switch(boost::math::fpclassify(b))
+   switch((::boost::math::fpclassify)(b))
    {
    case FP_NAN:
    case FP_INFINITE:
@@ -1038,7 +1057,7 @@ inline typename boost::enable_if_c<boost::is_arithmetic<A>::value>::type eval_fd
       result = zero;
       return;
    }
-   switch(boost::math::fpclassify(a))
+   switch((::boost::math::fpclassify)(a))
    {
    case FP_NAN:
       result = zero;
@@ -1055,12 +1074,6 @@ inline typename boost::enable_if_c<boost::is_arithmetic<A>::value>::type eval_fd
       result = zero;
 }
 
-template <class T1, class T2, class T3, class T4>
-inline void eval_fma(T1& result, const T2& a, const T3& b, const T4& c)
-{
-}
-
-
 template <class T>
 inline void eval_trunc(T& result, const T& a)
 {
@@ -1512,11 +1525,11 @@ namespace multiprecision{
    using boost::math::isnormal;
 
    typedef ::boost::math::policies::policy<
-      ::boost::math::policies::domain_error<::boost::math::policies::errno_on_error>,
-      ::boost::math::policies::pole_error<::boost::math::policies::errno_on_error>,
-      ::boost::math::policies::overflow_error<::boost::math::policies::errno_on_error>,
-      ::boost::math::policies::evaluation_error<::boost::math::policies::errno_on_error>,
-      ::boost::math::policies::rounding_error<::boost::math::policies::errno_on_error>
+      ::boost::math::policies::domain_error< ::boost::math::policies::errno_on_error>,
+      ::boost::math::policies::pole_error< ::boost::math::policies::errno_on_error>,
+      ::boost::math::policies::overflow_error< ::boost::math::policies::errno_on_error>,
+      ::boost::math::policies::evaluation_error< ::boost::math::policies::errno_on_error>,
+      ::boost::math::policies::rounding_error< ::boost::math::policies::errno_on_error>
    > c99_error_policy;
 
    template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
@@ -1955,6 +1968,200 @@ inline typename enable_if_c<number_category<B>::value == number_kind_integer, nu
    eval_integer_sqrt(s.backend(), r.backend(), x.backend());
    return s;
 }
+//
+// fma:
+//
+
+namespace default_ops {
+
+   struct fma_func
+   {
+      template <class B, class T, class U, class V>
+      void operator()(B& result, const T& a, const U& b, const V& c)const
+      {
+         eval_multiply_add(result, a, b, c);
+      }
+   };
+
+
+}
+
+template <class Backend, class U, class V>
+inline typename enable_if<
+   mpl::and_<
+      mpl::bool_<number_category<number<Backend, et_on> >::value == number_kind_floating_point>,
+      mpl::or_<
+         is_number<U>,
+         is_number_expression<U>,
+         is_arithmetic<U>
+      >,
+      mpl::or_<
+         is_number<V>,
+         is_number_expression<V>,
+         is_arithmetic<V>
+      >
+   >,
+   detail::expression<detail::function, default_ops::fma_func, number<Backend, et_on>, U, V>
+>::type
+fma(const number<Backend, et_on>& a, const U& b, const V& c)
+{
+   return detail::expression<detail::function, default_ops::fma_func, number<Backend, et_on>, U, V>(
+      default_ops::fma_func(), a, b, c);
+}
+
+template <class tag, class Arg1, class Arg2, class Arg3, class Arg4, class U, class V>
+inline typename enable_if<
+   mpl::and_<
+   mpl::bool_<number_category<typename detail::expression<tag, Arg1, Arg2, Arg3, Arg4>::result_type >::value == number_kind_floating_point>,
+   mpl::or_<
+   is_number<U>,
+   is_number_expression<U>,
+   is_arithmetic<U>
+   >,
+   mpl::or_<
+   is_number<V>,
+   is_number_expression<V>,
+   is_arithmetic<V>
+   >
+   >,
+   detail::expression<detail::function, default_ops::fma_func, detail::expression<tag, Arg1, Arg2, Arg3, Arg4>, U, V>
+>::type
+fma(const detail::expression<tag, Arg1, Arg2, Arg3, Arg4>& a, const U& b, const V& c)
+{
+   return detail::expression<detail::function, default_ops::fma_func, detail::expression<tag, Arg1, Arg2, Arg3, Arg4>, U, V>(
+      default_ops::fma_func(), a, b, c);
+}
+
+template <class Backend, class U, class V>
+inline typename enable_if<
+   mpl::and_<
+   mpl::bool_<number_category<number<Backend, et_off> >::value == number_kind_floating_point>,
+   mpl::or_<
+   is_number<U>,
+   is_number_expression<U>,
+   is_arithmetic<U>
+   >,
+   mpl::or_<
+   is_number<V>,
+   is_number_expression<V>,
+   is_arithmetic<V>
+   >
+   >,
+   number<Backend, et_off>
+>::type
+fma(const number<Backend, et_off>& a, const U& b, const V& c)
+{
+   using default_ops::eval_multiply_add;
+   number<Backend, et_off> result;
+   eval_multiply_add(result.backend(), number<Backend, et_off>::canonical_value(a), number<Backend, et_off>::canonical_value(b), number<Backend, et_off>::canonical_value(c));
+   return BOOST_MP_MOVE(result);
+}
+
+template <class U, class Backend, class V>
+inline typename enable_if<
+   mpl::and_<
+      mpl::bool_<number_category<number<Backend, et_on> >::value == number_kind_floating_point>,
+      is_arithmetic<U>,
+      mpl::or_<
+         is_number<V>,
+         is_number_expression<V>,
+         is_arithmetic<V>
+      >
+   >,
+   detail::expression<detail::function, default_ops::fma_func, U, number<Backend, et_on>, V>
+>::type
+fma(const U& a, const number<Backend, et_on>& b, const V& c)
+{
+   return detail::expression<detail::function, default_ops::fma_func, U, number<Backend, et_on>, V>(
+      default_ops::fma_func(), a, b, c);
+}
+
+template <class U, class tag, class Arg1, class Arg2, class Arg3, class Arg4, class V>
+inline typename enable_if<
+   mpl::and_<
+      mpl::bool_<number_category<typename detail::expression<tag, Arg1, Arg2, Arg3, Arg4>::result_type >::value == number_kind_floating_point>,
+      is_arithmetic<U>,
+      mpl::or_<
+         is_number<V>,
+         is_number_expression<V>,
+         is_arithmetic<V>
+      >
+   >,
+   detail::expression<detail::function, default_ops::fma_func, U, detail::expression<tag, Arg1, Arg2, Arg3, Arg4>, V>
+>::type
+fma(const U& a, const detail::expression<tag, Arg1, Arg2, Arg3, Arg4>& b, const V& c)
+{
+   return detail::expression<detail::function, default_ops::fma_func, U, detail::expression<tag, Arg1, Arg2, Arg3, Arg4>, V>(
+      default_ops::fma_func(), a, b, c);
+}
+
+template <class U, class Backend, class V>
+inline typename enable_if<
+   mpl::and_<
+      mpl::bool_<number_category<number<Backend, et_off> >::value == number_kind_floating_point>,
+      is_arithmetic<U>,
+      mpl::or_<
+         is_number<V>,
+         is_number_expression<V>,
+         is_arithmetic<V>
+      >
+   >,
+   number<Backend, et_off>
+>::type
+fma(const U& a, const number<Backend, et_off>& b, const V& c)
+{
+   using default_ops::eval_multiply_add;
+   number<Backend, et_off> result;
+   eval_multiply_add(result.backend(), number<Backend, et_off>::canonical_value(a), number<Backend, et_off>::canonical_value(b), number<Backend, et_off>::canonical_value(c));
+   return BOOST_MP_MOVE(result);
+}
+
+template <class U, class V, class Backend>
+inline typename enable_if<
+   mpl::and_<
+   mpl::bool_<number_category<number<Backend, et_on> >::value == number_kind_floating_point>,
+      is_arithmetic<U>,
+      is_arithmetic<V>
+   >,
+   detail::expression<detail::function, default_ops::fma_func, U, V, number<Backend, et_on> >
+>::type
+fma(const U& a, const V& b, const number<Backend, et_on>& c)
+{
+   return detail::expression<detail::function, default_ops::fma_func, U, V, number<Backend, et_on> >(
+      default_ops::fma_func(), a, b, c);
+}
+
+template <class U, class V, class tag, class Arg1, class Arg2, class Arg3, class Arg4>
+inline typename enable_if<
+   mpl::and_<
+   mpl::bool_<number_category<typename detail::expression<tag, Arg1, Arg2, Arg3, Arg4>::result_type >::value == number_kind_floating_point>,
+      is_arithmetic<U>,
+      is_arithmetic<V>
+   >,
+   detail::expression<detail::function, default_ops::fma_func, U, V, detail::expression<tag, Arg1, Arg2, Arg3, Arg4> >
+>::type
+fma(const U& a, const V& b, const detail::expression<tag, Arg1, Arg2, Arg3, Arg4>& c)
+{
+   return detail::expression<detail::function, default_ops::fma_func, U, V, detail::expression<tag, Arg1, Arg2, Arg3, Arg4> >(
+      default_ops::fma_func(), a, b, c);
+}
+
+template <class U, class V, class Backend>
+inline typename enable_if<
+   mpl::and_<
+   mpl::bool_<number_category<number<Backend, et_off> >::value == number_kind_floating_point>,
+      is_arithmetic<U>,
+      is_arithmetic<V>
+   >,
+   number<Backend, et_off>
+>::type
+fma(const U& a, const V& b, const number<Backend, et_off>& c)
+{
+   using default_ops::eval_multiply_add;
+   number<Backend, et_off> result;
+   eval_multiply_add(result.backend(), number<Backend, et_off>::canonical_value(a), number<Backend, et_off>::canonical_value(b), number<Backend, et_off>::canonical_value(c));
+   return BOOST_MP_MOVE(result);
+}
 
 template <class B, expression_template_option ExpressionTemplates>
 inline typename enable_if_c<number_category<B>::value == number_kind_integer, number<B, ExpressionTemplates> >::type

include/boost/multiprecision/detail/integer_ops.hpp

@@ -469,7 +469,19 @@ inline typename enable_if<
          is_integral<V>
       >
    >,
-   detail::expression<detail::function, default_ops::powm_func, T, U, V> >::type 
+   typename mpl::if_<
+      is_no_et_number<T>, 
+      T,
+      typename mpl::if_<
+         is_no_et_number<U>,
+         U,
+         typename mpl::if_<
+            is_no_et_number<V>,
+            V,
+            detail::expression<detail::function, default_ops::powm_func, T, U, V> >::type
+         >::type
+      >::type
+   >::type
    powm(const T& b, const U& p, const V& mod)
 {
    return detail::expression<detail::function, default_ops::powm_func, T, U, V>(

include/boost/multiprecision/detail/number_base.hpp

@@ -74,6 +74,18 @@ struct is_number : public mpl::false_ {};
 template <class Backend, expression_template_option ExpressionTemplates>
 struct is_number<number<Backend, ExpressionTemplates> > : public mpl::true_ {};
 
+template <class T>
+struct is_et_number : public mpl::false_ {};
+
+template <class Backend>
+struct is_et_number<number<Backend, et_on> > : public mpl::true_ {};
+
+template <class T>
+struct is_no_et_number : public mpl::false_ {};
+
+template <class Backend>
+struct is_no_et_number<number<Backend, et_off> > : public mpl::true_ {};
+
 namespace detail{
 
 // Forward-declare an expression wrapper
@@ -696,11 +708,11 @@ struct expression
    typedef typename right_middle_type::result_type right_middle_result_type;
    typedef typename right_type::result_type right_result_type;
    typedef typename combine_expression<
+      left_result_type,
       typename combine_expression<
-         typename combine_expression<left_result_type, left_middle_result_type>::type,
-         right_middle_result_type
-      >::type,
-      right_result_type
+         left_middle_result_type,
+         typename combine_expression<right_middle_result_type, right_result_type>::type
+      >::type
    >::type result_type;
    typedef tag tag_type;
 

include/boost/multiprecision/float128.hpp

@@ -489,6 +489,10 @@ inline void eval_atan2(float128_backend& result, const float128_backend& a, cons
 {
    result.value() = atan2q(a.value(), b.value());
 }
+inline void eval_multiply_add(float128_backend& result, const float128_backend& a, const float128_backend& b, const float128_backend& c)
+{
+   result.value() = fmaq(a.value(), b.value(), c.value());
+}
 
 inline std::size_t hash_value(const float128_backend& val)
 {

include/boost/multiprecision/mpfr.hpp

@@ -1452,6 +1452,31 @@ inline void eval_fmod(mpfr_float_backend<Digits10, AllocateType>& result, const
    mpfr_fmod(result.data(), a.data(), b.data(), GMP_RNDN);
 }
 
+template <unsigned Digits10, mpfr_allocation_type AllocateType>
+inline void eval_multiply_add(mpfr_float_backend<Digits10, AllocateType>& result, const mpfr_float_backend<Digits10, AllocateType>& a, const mpfr_float_backend<Digits10, AllocateType>& b)
+{
+   mpfr_fma(result.data(), a.data(), b.data(), result.data(), GMP_RNDN);
+}
+
+template <unsigned Digits10, mpfr_allocation_type AllocateType>
+inline void eval_multiply_add(mpfr_float_backend<Digits10, AllocateType>& result, const mpfr_float_backend<Digits10, AllocateType>& a, const mpfr_float_backend<Digits10, AllocateType>& b, const mpfr_float_backend<Digits10, AllocateType>& c)
+{
+   mpfr_fma(result.data(), a.data(), b.data(), c.data(), GMP_RNDN);
+}
+
+template <unsigned Digits10, mpfr_allocation_type AllocateType>
+inline void eval_multiply_subtract(mpfr_float_backend<Digits10, AllocateType>& result, const mpfr_float_backend<Digits10, AllocateType>& a, const mpfr_float_backend<Digits10, AllocateType>& b)
+{
+   mpfr_fms(result.data(), a.data(), b.data(), result.data(), GMP_RNDN);
+   result.negate();
+}
+
+template <unsigned Digits10, mpfr_allocation_type AllocateType>
+inline void eval_multiply_subtract(mpfr_float_backend<Digits10, AllocateType>& result, const mpfr_float_backend<Digits10, AllocateType>& a, const mpfr_float_backend<Digits10, AllocateType>& b, const mpfr_float_backend<Digits10, AllocateType>& c)
+{
+   mpfr_fms(result.data(), a.data(), b.data(), c.data(), GMP_RNDN);
+}
+
 template <unsigned Digits10, mpfr_allocation_type AllocateType>
 inline std::size_t hash_value(const mpfr_float_backend<Digits10, AllocateType>& val)
 {

test/test_arithmetic.hpp

@@ -742,7 +742,7 @@ void test_float_funcs(const boost::mpl::true_&)
    //
    // Test variable reuse in function calls, see https://svn.boost.org/trac/boost/ticket/8326
    //
-   Real a(2), b(10), c;
+   Real a(2), b(10), c, d;
    a = pow(a, b);
    BOOST_CHECK_EQUAL(a, 1024);
    a = 2;
@@ -898,6 +898,25 @@ void test_float_funcs(const boost::mpl::true_&)
    a = 4;
    b = atan2(a, b);
    BOOST_CHECK_CLOSE_FRACTION(b, Real(atan2(Real(4), Real(2))), tol);
+
+   // fma:
+   a = 2;
+   b = 4;
+   c = 6;
+   BOOST_CHECK_EQUAL(fma(a, b, c), 14);
+   BOOST_CHECK_EQUAL(fma(a, 4, c), 14);
+   BOOST_CHECK_EQUAL(fma(a, b, 6), 14);
+   BOOST_CHECK_EQUAL(fma(a, 4, 6), 14);
+   BOOST_CHECK_EQUAL(fma(a + 0, b, c), 14);
+   BOOST_CHECK_EQUAL(fma(a - 0, 4, c), 14);
+   BOOST_CHECK_EQUAL(fma(a * 1, b, 6), 14);
+   BOOST_CHECK_EQUAL(fma(a / 1, 4, 6), 14);
+   BOOST_CHECK_EQUAL(fma(2, b, c), 14);
+   BOOST_CHECK_EQUAL(fma(2, b, 6), 14);
+   BOOST_CHECK_EQUAL(fma(2, b * 1, c), 14);
+   BOOST_CHECK_EQUAL(fma(2, b + 0, 6), 14);
+   BOOST_CHECK_EQUAL(fma(2, 4, c), 14);
+   BOOST_CHECK_EQUAL(fma(2, 4, c + 0), 14);
 }
 
 template <class T, class U>