diff --git a/stan/math/prim.hpp b/stan/math/prim.hpp
index 208818ae116..1cbd9e8045a 100644
--- a/stan/math/prim.hpp
+++ b/stan/math/prim.hpp
@@ -14,6 +14,8 @@
 #include <stan/math/prim/err.hpp>
 
 #include <stan/math/prim/fun.hpp>
+#include <stan/math/prim/fun/is_even.hpp>
+#include <stan/math/prim/fun/is_odd.hpp>
 #include <stan/math/prim/functor.hpp>
 #include <stan/math/prim/prob.hpp>
 
diff --git a/stan/math/prim/fun/is_even.hpp b/stan/math/prim/fun/is_even.hpp
new file mode 100644
index 00000000000..2a50bfff887
--- /dev/null
+++ b/stan/math/prim/fun/is_even.hpp
@@ -0,0 +1,47 @@
+#ifndef STAN_MATH_PRIM_FUN_IS_EVEN_HPP
+#define STAN_MATH_PRIM_FUN_IS_EVEN_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/functor/apply_scalar_unary.hpp>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return 1 if the integer is even and 0 otherwise.
+ *
+ * @tparam T Integer type
+ * @param n Integer to check
+ * @return 1 if n is even, 0 otherwise
+ */
+template <typename T, require_integral_t<T>* = nullptr>
+inline int is_even(T n) {
+  return n % 2 == 0 ? 1 : 0;
+}
+
+/**
+ * Vectorization wrapper for is_even.
+ */
+struct is_even_fun {
+  template <typename T>
+  static inline int fun(const T& x) {
+    return is_even(static_cast<int>(x));
+  }
+};
+
+/**
+ * Apply is_even elementwise to a container.
+ *
+ * @tparam T Container type
+ * @param x Container of integers
+ * @return Elementwise parity results
+ */
+template <typename T, require_container_t<T>* = nullptr>
+inline auto is_even(T&& x) {
+  return apply_scalar_unary<is_even_fun, T>::apply(std::forward<T>(x));
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/prim/fun/is_odd.hpp b/stan/math/prim/fun/is_odd.hpp
new file mode 100644
index 00000000000..191f3f3cd59
--- /dev/null
+++ b/stan/math/prim/fun/is_odd.hpp
@@ -0,0 +1,47 @@
+#ifndef STAN_MATH_PRIM_FUN_IS_ODD_HPP
+#define STAN_MATH_PRIM_FUN_IS_ODD_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/functor/apply_scalar_unary.hpp>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return 1 if the integer is odd and 0 otherwise.
+ *
+ * @tparam T Integer type
+ * @param n Integer to check
+ * @return 1 if n is odd, 0 otherwise
+ */
+template <typename T, require_integral_t<T>* = nullptr>
+inline int is_odd(T n) {
+  return n % 2 != 0 ? 1 : 0;
+}
+
+/**
+ * Vectorization wrapper for is_odd.
+ */
+struct is_odd_fun {
+  template <typename T>
+  static inline int fun(const T& x) {
+    return is_odd(static_cast<int>(x));
+  }
+};
+
+/**
+ * Apply is_odd elementwise to a container.
+ *
+ * @tparam T Container type
+ * @param x Container of integers
+ * @return Elementwise parity results
+ */
+template <typename T, require_container_t<T>* = nullptr>
+inline auto is_odd(T&& x) {
+  return apply_scalar_unary<is_odd_fun, T>::apply(std::forward<T>(x));
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/test/unit/math/prim/fun/is_even_odd_test.cpp b/test/unit/math/prim/fun/is_even_odd_test.cpp
new file mode 100644
index 00000000000..f9c3f0657cb
--- /dev/null
+++ b/test/unit/math/prim/fun/is_even_odd_test.cpp
@@ -0,0 +1,69 @@
+#include <stan/math/prim.hpp>
+#include <gtest/gtest.h>
+#include <vector>
+
+TEST(MathFunctions, is_even_scalar) {
+  EXPECT_EQ(1, stan::math::is_even(0));
+  EXPECT_EQ(0, stan::math::is_even(1));
+  EXPECT_EQ(1, stan::math::is_even(2));
+  EXPECT_EQ(0, stan::math::is_even(3));
+
+  EXPECT_EQ(0, stan::math::is_even(-1));
+  EXPECT_EQ(1, stan::math::is_even(-2));
+  EXPECT_EQ(0, stan::math::is_even(-3));
+}
+
+TEST(MathFunctions, is_odd_scalar) {
+  EXPECT_EQ(0, stan::math::is_odd(0));
+  EXPECT_EQ(1, stan::math::is_odd(1));
+  EXPECT_EQ(0, stan::math::is_odd(2));
+  EXPECT_EQ(1, stan::math::is_odd(3));
+
+  EXPECT_EQ(1, stan::math::is_odd(-1));
+  EXPECT_EQ(0, stan::math::is_odd(-2));
+  EXPECT_EQ(1, stan::math::is_odd(-3));
+}
+
+TEST(MathFunctions, is_even_vector) {
+  std::vector<int> input{0, 1, 2, 3, -3, -4};
+  std::vector<int> expected{1, 0, 1, 0, 0, 1};
+
+  EXPECT_EQ(expected, stan::math::is_even(input));
+}
+
+TEST(MathFunctions, is_odd_vector) {
+  std::vector<int> input{0, 1, 2, 3, -3, -4};
+  std::vector<int> expected{0, 1, 0, 1, 1, 0};
+
+  EXPECT_EQ(expected, stan::math::is_odd(input));
+}
+
+TEST(MathFunctions, is_even_nested_vector) {
+  std::vector<std::vector<int>> input{{1, 2}, {3, 4}};
+
+  std::vector<std::vector<int>> expected{{0, 1}, {0, 1}};
+
+  EXPECT_EQ(expected, stan::math::is_even(input));
+}
+
+TEST(MathFunctions, is_odd_nested_vector) {
+  std::vector<std::vector<int>> input{{1, 2}, {3, 4}};
+
+  std::vector<std::vector<int>> expected{{1, 0}, {1, 0}};
+
+  EXPECT_EQ(expected, stan::math::is_odd(input));
+}
+
+TEST(MathFunctions, is_even_empty_vector) {
+  std::vector<int> input;
+  std::vector<int> expected;
+
+  EXPECT_EQ(expected, stan::math::is_even(input));
+}
+
+TEST(MathFunctions, is_odd_empty_vector) {
+  std::vector<int> input;
+  std::vector<int> expected;
+
+  EXPECT_EQ(expected, stan::math::is_odd(input));
+}