diff --git a/stan/math/prim/fun/square.hpp b/stan/math/prim/fun/square.hpp
index 0ed671bb96d..292f062a669 100644
--- a/stan/math/prim/fun/square.hpp
+++ b/stan/math/prim/fun/square.hpp
@@ -25,7 +25,12 @@ namespace math {
  */
 template <typename T, require_arithmetic_t<T>* = nullptr>
 inline double square(const T x) {
-  return std::pow(x, 2);
+  // widen to double first, then multiply: identical to the previous
+  // std::pow(x, 2) (a correctly-rounded square equals the rounded
+  // product), including the promoted-to-double semantics for integral
+  // arguments (where x * x could overflow), without the libm call.
+  const double x_d = x;
+  return x_d * x_d;
 }
 
 /**
diff --git a/stan/math/rev/fun/squared_distance.hpp b/stan/math/rev/fun/squared_distance.hpp
index 584114ee2c7..5ad19d9a8d5 100644
--- a/stan/math/rev/fun/squared_distance.hpp
+++ b/stan/math/rev/fun/squared_distance.hpp
@@ -21,11 +21,11 @@ namespace math {
 inline var squared_distance(const var& a, const var& b) {
   check_finite("squared_distance", "a", a);
   check_finite("squared_distance", "b", b);
-  return make_callback_vari(std::pow(a.val() - b.val(), 2),
-                            [a, b](const auto& vi) mutable {
-                              const double diff = 2.0 * (a.val() - b.val());
-                              a.adj() += vi.adj_ * diff;
-                              b.adj() -= vi.adj_ * diff;
+  const double diff = a.val() - b.val();
+  return make_callback_vari(diff * diff, [a, b, diff](const auto& vi) mutable {
+                              const double adj_diff = 2.0 * diff;
+                              a.adj() += vi.adj_ * adj_diff;
+                              b.adj() -= vi.adj_ * adj_diff;
                             });
 }
 
@@ -35,9 +35,9 @@ inline var squared_distance(const var& a, const var& b) {
 inline var squared_distance(const var& a, double b) {
   check_finite("squared_distance", "a", a);
   check_finite("squared_distance", "b", b);
-  return make_callback_vari(std::pow(a.val() - b, 2),
-                            [a, b](const auto& vi) mutable {
-                              a.adj() += vi.adj_ * 2.0 * (a.val() - b);
+  const double diff = a.val() - b;
+  return make_callback_vari(diff * diff, [a, b, diff](const auto& vi) mutable {
+                              a.adj() += vi.adj_ * 2.0 * diff;
                             });
 }