diff --git a/include/boost/math/special_functions/erf.hpp b/include/boost/math/special_functions/erf.hpp
index c1d9061db..dffb0c59c 100644
--- a/include/boost/math/special_functions/erf.hpp
+++ b/include/boost/math/special_functions/erf.hpp
@@ -153,7 +153,19 @@ T erf_imp(T z, bool invert, const Policy& pol, const Tag& t)
    if ((boost::math::isnan)(z))
       return policies::raise_domain_error("boost::math::erf<%1%>(%1%)", "Expected a finite argument but got %1%", z, pol);
 
-   if(z < 0)
+   if (fabs(z) < tools::root_epsilon<T>())
+   {
+      // Series[Erf[x], {x, 0, 4}]
+      // Series[Erfc[x], {x, 0, 4}]
+
+      const T term2 { z * 2 / constants::root_pi<T>() };
+
+      return invert ? 1 - term2 : term2;
+   }
+
+   const bool signbit_result = ((boost::math::signbit)(z) != 0);
+
+   if (signbit_result)
    {
       if(!invert)
          return -erf_imp(T(-z), invert, pol, t);
@@ -172,32 +184,32 @@ T erf_imp(T z, bool invert, const Policy& pol, const Tag& t)
    }
    else
    {
-      T x = z * z;
+      const T z_sq { z * z };
+
       if(z < 1.3f)
       {
          // Compute P:
          // This is actually good for z p to 2 or so, but the cutoff given seems
-         // to be the best compromise.  Performance wise, this is way quicker than anything else...
+         // to be the best compromise.  Regarding performance, this is way quicker than anything else...
          result = erf_series_near_zero_sum(z, pol);
       }
-      else if(x > 1 / tools::epsilon<T>())
+      else if(z_sq > 1 / tools::epsilon<T>())
       {
          // http://functions.wolfram.com/06.27.06.0006.02
          invert = !invert;
-         result = exp(-x) / (constants::root_pi<T>() * z);
+         result = exp(-z_sq) / (constants::root_pi<T>() * z);
       }
       else
       {
          // Compute Q:
          invert = !invert;
-         result = z * exp(-x);
+         result = z * exp(-z_sq);
          result /= boost::math::constants::root_pi<T>();
-         result *= upper_gamma_fraction(T(0.5f), x, policies::get_epsilon<T, Policy>());
+         result *= upper_gamma_fraction(T(0.5f), z_sq, policies::get_epsilon<T, Policy>());
       }
    }
-   if(invert)
-      result = 1 - result;
-   return result;
+
+   return ((!invert) ? result : 1 - result);
 }
 // LCOV_EXCL_STOP
 
diff --git a/test/test_erf.hpp b/test/test_erf.hpp
index b70c73953..2147938aa 100644
--- a/test/test_erf.hpp
+++ b/test/test_erf.hpp
@@ -174,6 +174,20 @@ void test_erf(T, const char* name)
       do_test_erfc_inv<T>(erfc_inv_big_data, name, "Inverse Erfc Function: extreme values");
    }
 
+   BOOST_CHECK_EQUAL(boost::math::erf(T(0)), T(0));
+   BOOST_CHECK_EQUAL(boost::math::erfc(T(0)), T(1));
+
+   BOOST_CHECK(boost::math::erf(boost::math::tools::root_epsilon<T>() / 8) > T(0));
+   BOOST_CHECK(boost::math::erf(boost::math::tools::epsilon<T>() / 32) >= T(0));
+   BOOST_CHECK(boost::math::erfc(boost::math::tools::epsilon<T>() / 32) <= T(1));
+
+   const bool has_negative_zero { ((boost::math::signbit)(-T(0)) != 0) };
+
+   if(has_negative_zero)
+   {
+      BOOST_CHECK_EQUAL(boost::math::erf(-T(0)), -T(0));
+   }
+
    BOOST_IF_CONSTEXPR(std::numeric_limits<T>::has_quiet_NaN)
    {
       BOOST_CHECK_THROW(boost::math::erf(std::numeric_limits<T>::quiet_NaN()), std::domain_error);