diff --git a/src/lib/math/pcurves/pcurves.cpp b/src/lib/math/pcurves/pcurves.cpp
index 9db4e189177..30b0dc06d7a 100644
--- a/src/lib/math/pcurves/pcurves.cpp
+++ b/src/lib/math/pcurves/pcurves.cpp
@@ -21,6 +21,13 @@ std::shared_ptr<const PrimeOrderCurve> PCurveInstance::secp192r1() {
 }
 #endif
 
+#if !defined(BOTAN_HAS_PCURVES_SECP224R1)
+//static
+std::shared_ptr<const PrimeOrderCurve> PCurveInstance::secp224r1() {
+   return nullptr;
+}
+#endif
+
 #if !defined(BOTAN_HAS_PCURVES_SECP256R1)
 //static
 std::shared_ptr<const PrimeOrderCurve> PCurveInstance::secp256r1() {
@@ -88,6 +95,8 @@ std::shared_ptr<const PrimeOrderCurve> PrimeOrderCurve::from_id(PrimeOrderCurveI
    switch(id.code()) {
       case PrimeOrderCurveId::secp192r1:
          return PCurveInstance::secp192r1();
+      case PrimeOrderCurveId::secp224r1:
+         return PCurveInstance::secp224r1();
       case PrimeOrderCurveId::secp256r1:
          return PCurveInstance::secp256r1();
       case PrimeOrderCurveId::secp384r1:
@@ -113,6 +122,7 @@ std::shared_ptr<const PrimeOrderCurve> PrimeOrderCurve::from_id(PrimeOrderCurveI
 std::vector<PrimeOrderCurveId> PrimeOrderCurveId::all() {
    return {
       PrimeOrderCurveId::secp192r1,
+      PrimeOrderCurveId::secp224r1,
       PrimeOrderCurveId::secp256r1,
       PrimeOrderCurveId::secp384r1,
       PrimeOrderCurveId::secp521r1,
@@ -129,6 +139,8 @@ std::string PrimeOrderCurveId::to_string() const {
    switch(this->code()) {
       case PrimeOrderCurveId::secp192r1:
          return "secp192r1";
+      case PrimeOrderCurveId::secp224r1:
+         return "secp224r1";
       case PrimeOrderCurveId::secp256r1:
          return "secp256r1";
       case PrimeOrderCurveId::secp384r1:
@@ -156,6 +168,8 @@ std::string PrimeOrderCurveId::to_string() const {
 std::optional<PrimeOrderCurveId> PrimeOrderCurveId::from_string(std::string_view name) {
    if(name == "secp192r1") {
       return PCurve::PrimeOrderCurveId::secp192r1;
+   } else if(name == "secp224r1") {
+      return PCurve::PrimeOrderCurveId::secp224r1;
    } else if(name == "secp256r1") {
       return PCurve::PrimeOrderCurveId::secp256r1;
    } else if(name == "secp384r1") {
diff --git a/src/lib/math/pcurves/pcurves_id.h b/src/lib/math/pcurves/pcurves_id.h
index 00fffb0b834..837f564357b 100644
--- a/src/lib/math/pcurves/pcurves_id.h
+++ b/src/lib/math/pcurves/pcurves_id.h
@@ -29,6 +29,8 @@ class BOTAN_TEST_API PrimeOrderCurveId final {
       enum class Code : uint8_t {
          /// secp192r1 aka P-192
          secp192r1,
+         /// secp224r1 aka P-224
+         secp224r1,
          /// secp256r1 aka P-256
          secp256r1,
          /// secp384r1 aka P-384
diff --git a/src/lib/math/pcurves/pcurves_impl/pcurves_impl.h b/src/lib/math/pcurves/pcurves_impl/pcurves_impl.h
index 33225014f85..5214fe108f8 100644
--- a/src/lib/math/pcurves/pcurves_impl/pcurves_impl.h
+++ b/src/lib/math/pcurves/pcurves_impl/pcurves_impl.h
@@ -72,7 +72,6 @@ class IntMod final {
       typedef typename Rep::W W;
 
       static constexpr auto P_MINUS_2 = p_minus<2>(P);
-      static constexpr auto P_PLUS_1_OVER_4 = p_plus_1_over_4(P);
 
    public:
       static constexpr size_t BITS = count_bits(P);
@@ -176,7 +175,7 @@ class IntMod final {
          return Self(Rep::redc(z));
       }
 
-      void square_n(size_t n) {
+      constexpr void square_n(size_t n) {
          std::array<W, 2 * N> z;
          for(size_t i = 0; i != n; ++i) {
             comba_sqr<N>(z.data(), this->data());
@@ -244,16 +243,54 @@ class IntMod final {
 
       /**
       * Return the modular square root, or zero if no root exists
-      *
-      * Current impl assumes p == 3 (mod 4)
       */
-      constexpr std::pair<Self, CT::Choice> sqrt() const
-         requires(Self::P_MOD_4 == 3)
-      {
-         auto z = pow_vartime(Self::P_PLUS_1_OVER_4);
-         const CT::Choice correct = (z.square() == *this);
-         z.conditional_assign(!correct, Self::zero());
-         return {z, correct};
+      constexpr std::pair<Self, CT::Choice> sqrt() const {
+         if constexpr(Self::P_MOD_4 == 3) {
+            constexpr auto P_PLUS_1_OVER_4 = p_plus_1_over_4(P);
+            auto z = pow_vartime(P_PLUS_1_OVER_4);
+            const CT::Choice correct = (z.square() == *this);
+            z.conditional_assign(!correct, Self::zero());
+            return {z, correct};
+         } else {
+            // Shanks-Tonelli, following I.4 in RFC 9380
+
+            /*
+            Constants:
+            1. c1, the largest integer such that 2^c1 divides q - 1.
+            2. c2 = (q - 1) / (2^c1)        # Integer arithmetic
+            3. c3 = (c2 - 1) / 2            # Integer arithmetic
+            4. c4, a non-square value in F
+            5. c5 = c4^c2 in F
+            */
+            constexpr size_t C1 = shanks_tonelli_c1(Self::P);
+            constexpr std::array<W, N> C2 = shanks_tonelli_c2(Self::P, C1);
+            constexpr std::array<W, N> C3 = shanks_tonelli_c3(C2);
+            constexpr std::array<W, N> P_MINUS_1_OVER_2 = p_minus_1_over_2(Self::P);
+            constexpr Self C4 = shanks_tonelli_c4<Self>(P_MINUS_1_OVER_2);
+            constexpr Self C5 = C4.pow_vartime(C2);
+
+            const Self& x = (*this);
+
+            auto z = x.pow_vartime(C3);
+            auto t = z.square();
+            t *= x;
+            z *= x;
+            auto b = t;
+            auto c = C5;
+
+            for(size_t i = C1; i >= 2; i--) {
+               b.square_n(i - 2);
+               const auto e = b.is_one();
+               z.conditional_assign(!e, z * c);
+               c.square_n(1);
+               t.conditional_assign(!e, t * c);
+               b = t;
+            }
+
+            const CT::Choice correct = (z.square() == *this);
+            z.conditional_assign(!correct, Self::zero());
+            return {z, correct};
+         }
       }
 
       constexpr CT::Choice operator==(const Self& other) const {
@@ -301,7 +338,7 @@ class IntMod final {
       }
 
       // Returns nullopt if the input is an encoding greater than or equal P
-      constexpr static std::optional<Self> deserialize(std::span<const uint8_t> bytes) {
+      static std::optional<Self> deserialize(std::span<const uint8_t> bytes) {
          // We could allow either short inputs or longer zero padded
          // inputs here, however it seems best to avoid non-canonical
          // representations unless required
@@ -474,7 +511,7 @@ class AffineCurvePoint {
 
       static constexpr FieldElement x3_ax_b(const FieldElement& x) { return (x.square() + Self::A) * x + Self::B; }
 
-      static constexpr std::optional<Self> deserialize(std::span<const uint8_t> bytes) {
+      static std::optional<Self> deserialize(std::span<const uint8_t> bytes) {
          if(bytes.size() == Self::BYTES) {
             if(bytes[0] != 0x04) {
                return {};
diff --git a/src/lib/math/pcurves/pcurves_impl/pcurves_util.h b/src/lib/math/pcurves/pcurves_impl/pcurves_util.h
index 21f4a627dc4..aa4c537166e 100644
--- a/src/lib/math/pcurves/pcurves_impl/pcurves_util.h
+++ b/src/lib/math/pcurves/pcurves_impl/pcurves_util.h
@@ -181,6 +181,66 @@ inline consteval std::array<W, N> p_minus_1_over_2(const std::array<W, N>& p) {
    return r;
 }
 
+template <WordType W, size_t N>
+inline consteval size_t shanks_tonelli_c1(const std::array<W, N>& p) {
+   auto get_bit = [&](size_t i) {
+      const size_t w = i / WordInfo<W>::bits;
+      const size_t b = i % WordInfo<W>::bits;
+      return static_cast<uint8_t>((p[w] >> b) & 0x01);
+   };
+
+   // Assumes low bit of P is set!
+   size_t b = 1;
+
+   while(get_bit(b) == 0) {
+      b += 1;
+   }
+
+   return b;
+}
+
+template <WordType W, size_t N>
+inline consteval std::array<W, N> shanks_tonelli_c2(const std::array<W, N>& p, size_t c1) {
+   auto c2 = p;
+   while(c1 > 0) {
+      shift_right<1>(c2);
+      c1--;
+   }
+   std::reverse(c2.begin(), c2.end());
+
+   return c2;
+}
+
+template <WordType W, size_t N>
+inline consteval std::array<W, N> shanks_tonelli_c3(const std::array<W, N>& c2) {
+   auto c3 = c2;
+   std::reverse(c3.begin(), c3.end());
+   shift_right<1>(c3);
+   std::reverse(c3.begin(), c3.end());
+   return c3;
+}
+
+template <typename Z, WordType W, size_t N>
+consteval auto shanks_tonelli_c4(const std::array<W, N>& p_minus_1_over_2) -> Z {
+   const auto one = Z::one();
+
+   // This is a silly performance hack; the first non-quadratic root in P-224
+   // is 11 so if we start the search there we save a little time.
+   auto z = Z::from_word(11);
+
+   for(;;) {
+      auto c = z.pow_vartime(p_minus_1_over_2);
+
+      auto is_square = c.is_zero() || c.is_one();
+
+      if(!is_square.as_bool()) {
+         return z;
+      }
+
+      z = z + one;
+   }
+}
+
 template <WordType W, size_t N>
 inline consteval size_t count_bits(const std::array<W, N>& p) {
    auto get_bit = [&](size_t i) {
diff --git a/src/lib/math/pcurves/pcurves_instance.h b/src/lib/math/pcurves/pcurves_instance.h
index 1ac05841a13..654e00edd0a 100644
--- a/src/lib/math/pcurves/pcurves_instance.h
+++ b/src/lib/math/pcurves/pcurves_instance.h
@@ -23,6 +23,8 @@ class PCurveInstance final {
 
       static std::shared_ptr<const PrimeOrderCurve> secp192r1();
 
+      static std::shared_ptr<const PrimeOrderCurve> secp224r1();
+
       static std::shared_ptr<const PrimeOrderCurve> secp256r1();
 
       static std::shared_ptr<const PrimeOrderCurve> secp384r1();
diff --git a/src/lib/math/pcurves/pcurves_secp224r1/info.txt b/src/lib/math/pcurves/pcurves_secp224r1/info.txt
new file mode 100644
index 00000000000..0186e1b8d41
--- /dev/null
+++ b/src/lib/math/pcurves/pcurves_secp224r1/info.txt
@@ -0,0 +1,13 @@
+<defines>
+PCURVES_SECP224R1 -> 20240716
+</defines>
+
+<module_info>
+name -> "PCurve secp224r1"
+brief -> "secp224r1"
+type -> "Internal"
+</module_info>
+
+<requires>
+pcurves_impl
+</requires>
diff --git a/src/lib/math/pcurves/pcurves_secp224r1/pcurves_secp224r1.cpp b/src/lib/math/pcurves/pcurves_secp224r1/pcurves_secp224r1.cpp
new file mode 100644
index 00000000000..1f244c29351
--- /dev/null
+++ b/src/lib/math/pcurves/pcurves_secp224r1/pcurves_secp224r1.cpp
@@ -0,0 +1,42 @@
+/*
+* (C) 2024 Jack Lloyd
+*
+* Botan is released under the Simplified BSD License (see license.txt)
+*/
+
+#include <botan/internal/pcurves_instance.h>
+
+#include <botan/internal/pcurves_solinas.h>
+#include <botan/internal/pcurves_wrap.h>
+
+namespace Botan::PCurve {
+
+namespace {
+
+namespace secp224r1 {
+
+// TODO Secp224r1Rep
+
+// clang-format off
+class Params final : public EllipticCurveParameters<
+   "FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF000000000000000000000001",
+   "FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFFFFFFFFFFFFFFFFFFFFE",
+   "B4050A850C04B3ABF54132565044B0B7D7BFD8BA270B39432355FFB4",
+   "FFFFFFFFFFFFFFFFFFFFFFFFFFFF16A2E0B8F03E13DD29455C5C2A3D",
+   "B70E0CBD6BB4BF7F321390B94A03C1D356C21122343280D6115C1D21",
+   "BD376388B5F723FB4C22DFE6CD4375A05A07476444D5819985007E34"> {
+};
+
+// clang-format on
+
+class Curve final : public EllipticCurve<Params> {};
+
+}  // namespace secp224r1
+
+}  // namespace
+
+std::shared_ptr<const PrimeOrderCurve> PCurveInstance::secp224r1() {
+   return PrimeOrderCurveImpl<secp224r1::Curve>::instance();
+}
+
+}  // namespace Botan::PCurve
diff --git a/src/scripts/ci_build.py b/src/scripts/ci_build.py
index 5171fbaa857..21e1bc90da0 100755
--- a/src/scripts/ci_build.py
+++ b/src/scripts/ci_build.py
@@ -252,7 +252,7 @@ def sanitize_kv(some_string):
                 'ed25519_sign', 'elgamal_decrypt', 'elgamal_encrypt', 'elgamal_keygen',
                 'ffi_dh', 'ffi_dsa', 'ffi_elgamal', 'frodo_kat_tests', 'hash_nist_mc',
                 'hss_lms_keygen', 'hss_lms_sign', 'mce_keygen', 'passhash9', 'pbkdf',
-                'pcurves_points', 'pwdhash', 'rsa_encrypt', 'rsa_pss', 'rsa_pss_raw', 'scrypt',
+                'pcurves_arith', 'pwdhash', 'rsa_encrypt', 'rsa_pss', 'rsa_pss_raw', 'scrypt',
                 'sphincsplus', 'sphincsplus_fors', 'sphincsplus_keygen', 'srp6_kat',
                 'srp6_rt', 'unit_tls', 'x509_path_bsi', 'x509_path_rsa_pss',
                 'xmss_keygen', 'xmss_keygen_reference', 'xmss_sign', 'xmss_unit_tests',
diff --git a/src/tests/test_pcurves.cpp b/src/tests/test_pcurves.cpp
index 16fb98123b4..7a5c1fcffdc 100644
--- a/src/tests/test_pcurves.cpp
+++ b/src/tests/test_pcurves.cpp
@@ -195,7 +195,7 @@ class Pcurve_Ecdsa_Sign_Tests final : public Text_Based_Test {
 
 BOTAN_REGISTER_TEST("pcurves", "pcurves_ecdsa_sign", Pcurve_Ecdsa_Sign_Tests);
 
-class Pcurve_Point_Tests final : public Test {
+class Pcurve_Arithmetic_Tests final : public Test {
    public:
       std::vector<Test::Result> run() override {
          std::vector<Test::Result> results;
@@ -339,7 +339,62 @@ class Pcurve_Point_Tests final : public Test {
       }
 };
 
-BOTAN_REGISTER_TEST("pcurves", "pcurves_points", Pcurve_Point_Tests);
+BOTAN_REGISTER_TEST("pcurves", "pcurves_arith", Pcurve_Arithmetic_Tests);
+
+class Pcurve_PointEnc_Tests final : public Test {
+   public:
+      std::vector<Test::Result> run() override {
+         std::vector<Test::Result> results;
+
+         auto& rng = Test::rng();
+
+         for(auto id : Botan::PCurve::PrimeOrderCurveId::all()) {
+            Test::Result result("Pcurves point operations " + id.to_string());
+
+            result.start_timer();
+
+            auto curve = Botan::PCurve::PrimeOrderCurve::from_id(id);
+
+            if(!curve) {
+               result.test_note("Skipping test due to missing pcurve " + id.to_string());
+               continue;
+            }
+
+            for(size_t trial = 0; trial != 100; ++trial) {
+               const auto scalar = curve->random_scalar(rng);
+               const auto pt = curve->mul_by_g(scalar, rng).to_affine();
+
+               const auto pt_u = pt.serialize();
+               result.test_eq("Expected uncompressed header", static_cast<size_t>(pt_u[0]), 0x04);
+               const size_t fe_bytes = (pt_u.size() - 1) / 2;
+               const auto pt_c = pt.serialize_compressed();
+
+               result.test_eq("Expected compressed size", pt_c.size(), 1 + fe_bytes);
+               result.confirm("Expected compressed header", pt_c[0] == 0x02 || pt_c[0] == 0x03);
+
+               if(auto d_pt_u = curve->deserialize_point(pt_u)) {
+                  result.test_eq("Deserializing uncompressed returned correct point", d_pt_u->serialize(), pt_u);
+               } else {
+                  result.test_failure("Failed to deserialize uncompressed point");
+               }
+
+               if(auto d_pt_c = curve->deserialize_point(pt_c)) {
+                  result.test_eq("Deserializing compressed returned correct point", d_pt_c->serialize(), pt_u);
+               } else {
+                  result.test_failure("Failed to deserialize compressed point");
+               }
+            }
+
+            result.end_timer();
+
+            results.push_back(result);
+         }
+
+         return results;
+      }
+};
+
+BOTAN_REGISTER_TEST("pcurves", "pcurves_point_enc", Pcurve_PointEnc_Tests);
 
 #endif