diff --git a/include/boost/multiprecision/mpfi.hpp b/include/boost/multiprecision/mpfi.hpp
index 537eb3ee..b2b370c6 100644
--- a/include/boost/multiprecision/mpfi.hpp
+++ b/include/boost/multiprecision/mpfi.hpp
@@ -477,7 +477,10 @@ inline void eval_subtract(mpfi_float_backend<D1>& result, const mpfi_float_backe
 template <unsigned D1, unsigned D2>
 inline void eval_multiply(mpfi_float_backend<D1>& result, const mpfi_float_backend<D2>& o)
 {
-   mpfi_mul(result.data(), result.data(), o.data());
+   if(&result == &o)
+      mpfi_sqr(result.data(), o.data());
+   else
+      mpfi_mul(result.data(), result.data(), o.data());
 }
 template <unsigned D1, unsigned D2>
 inline void eval_divide(mpfi_float_backend<D1>& result, const mpfi_float_backend<D2>& o)
@@ -609,7 +612,10 @@ inline void eval_subtract(mpfi_float_backend<D1>& a, long x, const mpfi_float_ba
 template <unsigned D1, unsigned D2, unsigned D3>
 inline void eval_multiply(mpfi_float_backend<D1>& a, const mpfi_float_backend<D2>& x, const mpfi_float_backend<D3>& y)
 {
-   mpfi_mul(a.data(), x.data(), y.data());
+   if(&x == &y)
+      mpfi_sqr(a.data(), x.data());
+   else
+      mpfi_mul(a.data(), x.data(), y.data());
 }
 template <unsigned D1, unsigned D2>
 inline void eval_multiply(mpfi_float_backend<D1>& a, const mpfi_float_backend<D2>& x, unsigned long y)
diff --git a/include/boost/multiprecision/mpfr.hpp b/include/boost/multiprecision/mpfr.hpp
index 18e7024b..da864fb4 100644
--- a/include/boost/multiprecision/mpfr.hpp
+++ b/include/boost/multiprecision/mpfr.hpp
@@ -952,7 +952,10 @@ inline void eval_subtract(mpfr_float_backend<D1, A1>& result, const mpfr_float_b
 template <unsigned D1, unsigned D2, mpfr_allocation_type A1, mpfr_allocation_type A2>
 inline void eval_multiply(mpfr_float_backend<D1, A1>& result, const mpfr_float_backend<D2, A2>& o)
 {
-   mpfr_mul(result.data(), result.data(), o.data(), GMP_RNDN);
+   if(&o == &result)
+      mpfr_sqr(result.data(), o.data(), GMP_RNDN);
+   else
+      mpfr_mul(result.data(), result.data(), o.data(), GMP_RNDN);
 }
 template <unsigned D1, unsigned D2, mpfr_allocation_type A1, mpfr_allocation_type A2>
 inline void eval_divide(mpfr_float_backend<D1, A1>& result, const mpfr_float_backend<D2, A2>& o)
@@ -1084,7 +1087,10 @@ inline void eval_subtract(mpfr_float_backend<D1, A1>& a, long x, const mpfr_floa
 template <unsigned D1, unsigned D2, mpfr_allocation_type A1, mpfr_allocation_type A2, unsigned D3>
 inline void eval_multiply(mpfr_float_backend<D1, A1>& a, const mpfr_float_backend<D2, A2>& x, const mpfr_float_backend<D3>& y)
 {
-   mpfr_mul(a.data(), x.data(), y.data(), GMP_RNDN);
+   if(&x == &y)
+      mpfr_sqr(a.data(), x.data(), GMP_RNDN);
+   else
+      mpfr_mul(a.data(), x.data(), y.data(), GMP_RNDN);
 }
 template <unsigned D1, unsigned D2, mpfr_allocation_type A1, mpfr_allocation_type A2>
 inline void eval_multiply(mpfr_float_backend<D1, A1>& a, const mpfr_float_backend<D2, A2>& x, unsigned long y)