diff --git a/include/xsimd/arch/common/xsimd_common_trigo.hpp b/include/xsimd/arch/common/xsimd_common_trigo.hpp
index 00ad75076..462c94026 100644
--- a/include/xsimd/arch/common/xsimd_common_trigo.hpp
+++ b/include/xsimd/arch/common/xsimd_common_trigo.hpp
@@ -358,11 +358,22 @@ namespace xsimd
             const batch_type q = abs(self / other);
             const batch_type q_p = abs(other / self);
 #else
-            const batch_type q = abs(self / other);
+            const auto both_zero = (self == batch_type(0.)) && (other == batch_type(0.));
+            auto use_unit_ratio = both_zero;
+#ifndef XSIMD_NO_INFINITIES
+            use_unit_ratio = use_unit_ratio || (isinf(self) && isinf(other));
+#endif
+            const batch_type q = abs(select(use_unit_ratio, batch_type(1.), self)
+                                     / select(use_unit_ratio, batch_type(1.), other));
             const batch_type q_p = 1. / q;
 #endif
             const batch_type z = detail::kernel_atan(q, q_p);
-            return select(other > batch_type(0.), z, constants::pi<batch_type>() - z) * signnz(self);
+            batch_type result = select(other > batch_type(0.), z, constants::pi<batch_type>() - z) * signnz(self);
+#ifndef __FAST_MATH__
+            const auto zero_angle = select(signnz(other) < batch_type(0.), constants::pi<batch_type>(), batch_type(0.));
+            result = select(both_zero, copysign(zero_angle, self), result);
+#endif
+            return result;
         }
 
         // cos
diff --git a/test/test_trigonometric.cpp b/test/test_trigonometric.cpp
index d3e8862a6..af1dc4d05 100644
--- a/test/test_trigonometric.cpp
+++ b/test/test_trigonometric.cpp
@@ -179,6 +179,52 @@ struct trigonometric_test
             }
         }
     }
+    void test_atan2_pairs()
+    {
+#ifndef __FAST_MATH__
+        // Fast-math omits zero and infinity recovery and returns NaN.
+        const value_type spec(0);
+        const value_type nspec = -spec;
+        const value_type big = std::numeric_limits<value_type>::infinity();
+#else
+        const value_type spec(0.5);
+        const value_type nspec = -spec;
+        const value_type big(8);
+#endif
+        const value_type ys[] = { -big, -value_type(3), -value_type(1), nspec, spec, value_type(1), value_type(3), big };
+        const value_type xs[] = { -big, -value_type(4), -value_type(1), nspec, spec, value_type(1), value_type(4), big };
+        vector_type y_values;
+        vector_type x_values;
+        y_values.reserve(64 + size);
+        x_values.reserve(64 + size);
+        for (value_type y : ys)
+        {
+            for (value_type x : xs)
+            {
+                y_values.push_back(y);
+                x_values.push_back(x);
+            }
+        }
+        while (y_values.size() % size != 0)
+        {
+            y_values.push_back(y_values.front());
+            x_values.push_back(x_values.front());
+        }
+        for (size_t i = 0; i < y_values.size(); i += size)
+        {
+            batch_type yb, xb, out, ref;
+            detail::load_batch(yb, y_values, i);
+            detail::load_batch(xb, x_values, i);
+            out = atan2(yb, xb);
+            vector_type scalar(size);
+            for (size_t lane = 0; lane < size; ++lane)
+            {
+                scalar[lane] = std::atan2(y_values[i + lane], x_values[i + lane]);
+            }
+            detail::load_batch(ref, scalar, 0);
+            CHECK_BATCH_EQ(ref, out);
+        }
+    }
 };
 
 TEST_CASE_TEMPLATE("[trigonometric]", B, BATCH_FLOAT_TYPES)
@@ -189,6 +235,11 @@ TEST_CASE_TEMPLATE("[trigonometric]", B, BATCH_FLOAT_TYPES)
         Test.test_trigonometric_functions();
     }
 
+    SUBCASE("atan2 pairs")
+    {
+        Test.test_atan2_pairs();
+    }
+
     SUBCASE("reciprocal")
     {
         Test.test_reciprocal_functions();