diff --git a/benches/const_monty.rs b/benches/const_monty.rs
index 136762c2..7405896b 100644
--- a/benches/const_monty.rs
+++ b/benches/const_monty.rs
@@ -193,10 +193,94 @@ fn bench_montgomery_ops<M: Measurement>(group: &mut BenchmarkGroup<'_, M>) {
     }
 }
 
+fn bench_montgomery_sqrt<M: Measurement>(group: &mut BenchmarkGroup<'_, M>) {
+    use crypto_bigint::{U256, const_prime_monty_params, modular::ConstPrimeMontyParams};
+
+    {
+        // P-256 field modulus
+        // p = 3 mod 4, s = 1, uses Shanks algorithm
+        const_prime_monty_params!(
+            P256Field,
+            U256,
+            "ffffffff00000001000000000000000000000000ffffffffffffffffffffffff"
+        );
+        assert_eq!(P256Field::PRIME_PARAMS.s().get(), 1);
+        type ConstForm = crypto_bigint::modular::ConstMontyForm<P256Field, { U256::LIMBS }>;
+
+        let mut rng = ChaCha8Rng::from_seed([7u8; 32]);
+        group.bench_function("sqrt, U256, s=1", |b| {
+            b.iter_batched(
+                || {
+                    let x =
+                        U256::random_mod_vartime(&mut rng, P256Field::PARAMS.modulus().as_nz_ref());
+                    ConstForm::new(&x)
+                },
+                |x| x.sqrt(),
+                BatchSize::SmallInput,
+            );
+        });
+    }
+
+    {
+        // P-256 scalar modulus
+        // p = 17 mod 32, s = 4
+        const_prime_monty_params!(
+            P256Scalar,
+            U256,
+            "ffffffff00000000ffffffffffffffffbce6faada7179e84f3b9cac2fc632551"
+        );
+        assert_eq!(P256Scalar::PRIME_PARAMS.s().get(), 4);
+        type ConstForm = crypto_bigint::modular::ConstMontyForm<P256Scalar, { U256::LIMBS }>;
+
+        let mut rng = ChaCha8Rng::from_seed([7u8; 32]);
+        group.bench_function("sqrt, U256, s=4", |b| {
+            b.iter_batched(
+                || {
+                    let x = U256::random_mod_vartime(
+                        &mut rng,
+                        P256Scalar::PARAMS.modulus().as_nz_ref(),
+                    );
+                    ConstForm::new(&x)
+                },
+                |x| x.sqrt(),
+                BatchSize::SmallInput,
+            );
+        });
+    }
+
+    {
+        // K-256 scalar modulus
+        // s = 6, uses Tonelli-Shanks
+        const_prime_monty_params!(
+            K256Scalar,
+            U256,
+            "FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEBAAEDCE6AF48A03BBFD25E8CD0364141"
+        );
+        assert_eq!(K256Scalar::PRIME_PARAMS.s().get(), 6);
+        type ConstForm = crypto_bigint::modular::ConstMontyForm<K256Scalar, { U256::LIMBS }>;
+
+        let mut rng = ChaCha8Rng::from_seed([7u8; 32]);
+        group.bench_function("sqrt, U256, s=6", |b| {
+            b.iter_batched(
+                || {
+                    let x = U256::random_mod_vartime(
+                        &mut rng,
+                        K256Scalar::PARAMS.modulus().as_nz_ref(),
+                    );
+                    ConstForm::new(&x)
+                },
+                |x| x.sqrt(),
+                BatchSize::SmallInput,
+            );
+        });
+    }
+}
+
 fn bench_montgomery(c: &mut Criterion) {
     let mut group = c.benchmark_group("Const Montgomery arithmetic");
     bench_montgomery_conversion(&mut group);
     bench_montgomery_ops(&mut group);
+    bench_montgomery_sqrt(&mut group);
     group.finish();
 }
 
diff --git a/src/modular.rs b/src/modular.rs
index 6434ed0e..8fb74ed6 100644
--- a/src/modular.rs
+++ b/src/modular.rs
@@ -28,16 +28,19 @@ mod div_by_2;
 mod monty_params;
 mod mul;
 mod pow;
+mod prime_params;
 pub(crate) mod safegcd;
+mod sqrt;
 mod sub;
 
 #[cfg(feature = "alloc")]
 pub(crate) mod boxed_monty_form;
 
 pub use self::{
-    const_monty_form::{ConstMontyForm, ConstMontyParams},
+    const_monty_form::{ConstMontyForm, ConstMontyParams, ConstPrimeMontyParams},
     fixed_monty_form::FixedMontyForm,
     monty_params::{FixedMontyParams, MontyParams},
+    prime_params::PrimeParams,
 };
 
 pub(crate) use self::safegcd::SafeGcdInverter;
diff --git a/src/modular/const_monty_form.rs b/src/modular/const_monty_form.rs
index e6844acd..bee62aad 100644
--- a/src/modular/const_monty_form.rs
+++ b/src/modular/const_monty_form.rs
@@ -8,10 +8,11 @@ mod mul;
 mod neg;
 mod pow;
 mod reduce;
+mod sqrt;
 mod sub;
 
 use super::{
-    FixedMontyParams, Retrieve, div_by_2::div_by_2, mul::mul_montgomery_form,
+    FixedMontyParams, PrimeParams, Retrieve, div_by_2::div_by_2, mul::mul_montgomery_form,
     reduction::montgomery_retrieve,
 };
 use crate::{ConstOne, ConstZero, CtEq, Odd, One, Uint, Zero};
@@ -47,6 +48,16 @@ pub trait ConstMontyParams<const LIMBS: usize>:
     const PARAMS: FixedMontyParams<LIMBS>;
 }
 
+/// Trait representing a prime modulus and its associated constants for converting in
+/// and out of Montgomery form.
+///
+/// To define a type which impls this trait, use the
+/// [`const_prime_monty_params!`][`crate::const_prime_monty_params`] macro.
+pub trait ConstPrimeMontyParams<const LIMBS: usize>: ConstMontyParams<LIMBS> {
+    /// Prime parameters constant.
+    const PRIME_PARAMS: PrimeParams<LIMBS>;
+}
+
 /// An integer in Montgomery form modulo `MOD`, represented using `LIMBS` limbs.
 /// The modulus is constant, so it cannot be set at runtime.
 ///
diff --git a/src/modular/const_monty_form/macros.rs b/src/modular/const_monty_form/macros.rs
index e7d4fc70..d62b9e4a 100644
--- a/src/modular/const_monty_form/macros.rs
+++ b/src/modular/const_monty_form/macros.rs
@@ -1,7 +1,7 @@
 //! [`ConstMontyForm`]/[`ConstMontyParams`] support macros.
 
 #[cfg(doc)]
-use crate::modular::{ConstMontyForm, ConstMontyParams};
+use crate::modular::{ConstMontyForm, ConstMontyParams, ConstPrimeMontyParams};
 
 /// Create a type representing a modulus which impls the [`ConstMontyParams`] trait with the given
 /// name, type, value (in big endian hex), and optional documentation string.
@@ -45,6 +45,52 @@ macro_rules! const_monty_params {
     };
 }
 
+/// Create a type representing a prime modulus which impls the [`ConstPrimeMontyParams`]
+/// trait with the given name, type, value (in big endian hex), and optional documentation
+/// string.
+///
+/// # Usage
+///
+/// ```
+/// use crypto_bigint::{U256, const_prime_monty_params};
+///
+/// const_prime_monty_params!(
+///     MyModulus,
+///     U256,
+///     "73eda753299d7d483339d80809a1d80553bda402fffe5bfeffffffff00000001",
+///     "Docs for my modulus"
+/// );
+/// ```
+///
+/// The modulus _must_ be odd and prime, or this will panic.
+#[macro_export]
+macro_rules! const_prime_monty_params {
+    ($name:ident, $uint_type:ty, $value:expr) => {
+        $crate::const_prime_monty_params!(
+            $name,
+            $uint_type,
+            $value,
+            "Modulus which impls `ConstPrimeMontyParams`"
+        );
+    };
+    ($name:ident, $uint_type:ty, $value:expr, $doc:expr) => {
+        $crate::const_monty_params!(
+            $name,
+            $uint_type,
+            $value,
+            "Modulus which impls `ConstPrimeMontyParams`"
+        );
+
+        impl $crate::modular::ConstPrimeMontyParams<{ <$uint_type>::LIMBS }> for $name {
+            const PRIME_PARAMS: $crate::modular::PrimeParams<{ <$uint_type>::LIMBS }> =
+                $crate::modular::PrimeParams::new_vartime(
+                    &<$name as $crate::modular::ConstMontyParams<{ <$uint_type>::LIMBS }>>::PARAMS,
+                )
+                .expect("cannot derive prime parameters");
+        }
+    };
+}
+
 /// Creates a type alias to [`ConstMontyForm`] with the given [`ConstMontyParams`].
 ///
 /// # Usage
diff --git a/src/modular/const_monty_form/sqrt.rs b/src/modular/const_monty_form/sqrt.rs
new file mode 100644
index 00000000..49f2116c
--- /dev/null
+++ b/src/modular/const_monty_form/sqrt.rs
@@ -0,0 +1,59 @@
+use core::marker::PhantomData;
+
+use super::ConstPrimeMontyParams;
+use crate::{
+    CtOption,
+    modular::{ConstMontyForm, sqrt::sqrt_montgomery_form},
+};
+
+impl<const LIMBS: usize, MOD> ConstMontyForm<MOD, LIMBS>
+where
+    MOD: ConstPrimeMontyParams<LIMBS>,
+{
+    /// Compute the modular square root for `self`, if it exists.
+    #[must_use]
+    pub const fn sqrt(&self) -> CtOption<Self> {
+        let res = sqrt_montgomery_form(self.as_montgomery(), &MOD::PARAMS, &MOD::PRIME_PARAMS);
+        let is_some = res.is_some();
+        CtOption::new(
+            Self {
+                montgomery_form: *res.as_inner_unchecked(),
+                phantom: PhantomData,
+            },
+            is_some,
+        )
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use crate::{
+        U256, const_prime_monty_params,
+        modular::{ConstMontyForm, ConstPrimeMontyParams},
+    };
+
+    #[test]
+    fn check_sqrt() {
+        // P-256 field modulus
+        // p = 3 mod 4, s = 1, uses Shanks algorithm
+        const_prime_monty_params!(
+            P256Field,
+            U256,
+            "ffffffff00000001000000000000000000000000ffffffffffffffffffffffff"
+        );
+        assert_eq!(P256Field::PRIME_PARAMS.s().get(), 1);
+        type ConstForm = ConstMontyForm<P256Field, { U256::LIMBS }>;
+
+        // four is square
+        let four_monty = ConstForm::new(&U256::from(4u32));
+        assert_eq!(
+            four_monty.sqrt().expect("ensured square"),
+            ConstForm::new(&U256::from(2u32))
+        );
+
+        // generator must be non-residue
+        let generator = U256::from_u32(P256Field::PRIME_PARAMS.generator().get());
+        let gen_monty = ConstForm::new(&generator);
+        assert!(gen_monty.sqrt().is_none().to_bool_vartime());
+    }
+}
diff --git a/src/modular/prime_params.rs b/src/modular/prime_params.rs
new file mode 100644
index 00000000..fc6b2d48
--- /dev/null
+++ b/src/modular/prime_params.rs
@@ -0,0 +1,246 @@
+//! Parameter calculation for prime moduli.
+
+use core::num::NonZeroU32;
+use ctutils::{CtAssignSlice, CtEqSlice, CtSelectUsingCtAssign};
+
+use super::{FixedMontyForm, FixedMontyParams};
+use crate::{Choice, CtAssign, CtEq, Odd, Uint};
+
+#[cfg(feature = "subtle")]
+use crate::CtSelect;
+
+/// Parameters for supporting efficient computations on integers in Montgomery form
+/// with a prime modulus.
+#[derive(Debug, Copy, Clone)]
+pub struct PrimeParams<const LIMBS: usize> {
+    /// A constant such that the modulus `p = t•2^s+1` for `s > 0` and some odd `t`.
+    pub(super) s: NonZeroU32,
+    /// The smallest primitive root of the modulus.
+    pub(super) generator: NonZeroU32,
+    /// The exponent to use in computing a modular square root.
+    pub(super) sqrt_exp: Uint<LIMBS>,
+    /// An s'th root of unity for the modulus, in Montgomery form.
+    pub(super) monty_root_unity: Uint<LIMBS>,
+    /// Equal to `monty_root_unity^2 mod p`.
+    pub(super) monty_root_unity_p2: Uint<LIMBS>,
+}
+
+impl<const LIMBS: usize> PrimeParams<LIMBS> {
+    /// Instantiates a new set of [`PrimeParams`] given [`FixedMontyParams`] for a prime modulus.
+    ///
+    /// This method will return `None` if the modulus is determined to be non-prime, however
+    /// this is not an exhaustive check and non-prime values can be accepted.
+    #[must_use]
+    #[allow(clippy::unwrap_in_result, clippy::missing_panics_doc)]
+    pub const fn new_vartime(params: &FixedMontyParams<LIMBS>) -> Option<Self> {
+        let p = params.modulus();
+        let p_minus_one = p.as_ref().set_bit_vartime(0, false);
+        let s = NonZeroU32::new(p_minus_one.trailing_zeros_vartime()).expect("ensured non-zero");
+
+        let Some((generator, gen_uint)) = find_primitive_root(p) else {
+            return None;
+        };
+
+        // if s=1 and p is a power of a prime then -1 is always a root of unity
+        let (exp, root) = if s.get() == 1 {
+            // (p+1)/4
+            let exp = p.as_ref().shr_vartime(2).wrapping_add(&Uint::ONE);
+            let root = FixedMontyForm::new(&p_minus_one, params);
+            (exp, root)
+        } else {
+            // t = (p-1)/2^s
+            let t = p.as_ref().shr_vartime(s.get());
+            // exp = (t-1)/2
+            let exp = t.shr_vartime(1);
+            // the s'th root of unity is calculated as `generator^t`
+            let root = FixedMontyForm::new(&gen_uint, params).pow_vartime(&t);
+            // root^(2^(s-1)) must be equal to -1
+            let check = root.square_repeat_vartime(s.get() - 1);
+            if !Uint::eq(&check.retrieve(), &p_minus_one).to_bool_vartime() {
+                return None;
+            }
+            (exp, root)
+        };
+
+        Some(Self {
+            s,
+            generator,
+            sqrt_exp: exp,
+            monty_root_unity: root.to_montgomery(),
+            monty_root_unity_p2: root.square().to_montgomery(),
+        })
+    }
+
+    /// Get the constant 'generator' used in modular square root calculation.
+    #[must_use]
+    pub const fn generator(&self) -> NonZeroU32 {
+        self.generator
+    }
+
+    /// Get the constant 's' used in modular square root calculation.
+    #[must_use]
+    pub const fn s(&self) -> NonZeroU32 {
+        self.s
+    }
+}
+
+impl<const LIMBS: usize> CtAssign for PrimeParams<LIMBS> {
+    fn ct_assign(&mut self, other: &Self, choice: Choice) {
+        self.s.ct_assign(&other.s, choice);
+        self.generator.ct_assign(&other.generator, choice);
+        self.sqrt_exp.ct_assign(&other.sqrt_exp, choice);
+        self.monty_root_unity
+            .ct_assign(&other.monty_root_unity, choice);
+        self.monty_root_unity_p2
+            .ct_assign(&other.monty_root_unity_p2, choice);
+    }
+}
+impl<const LIMBS: usize> CtAssignSlice for PrimeParams<LIMBS> {}
+impl<const LIMBS: usize> CtSelectUsingCtAssign for PrimeParams<LIMBS> {}
+
+impl<const LIMBS: usize> CtEq for PrimeParams<LIMBS> {
+    fn ct_eq(&self, other: &Self) -> Choice {
+        self.s.ct_eq(&other.s)
+            & self.generator.ct_eq(&other.generator)
+            & self.sqrt_exp.ct_eq(&other.sqrt_exp)
+            & self.monty_root_unity.ct_eq(&other.monty_root_unity)
+            & self.monty_root_unity_p2.ct_eq(&other.monty_root_unity_p2)
+    }
+}
+impl<const LIMBS: usize> CtEqSlice for PrimeParams<LIMBS> {}
+
+#[cfg(feature = "subtle")]
+impl<const LIMBS: usize> subtle::ConstantTimeEq for PrimeParams<LIMBS> {
+    fn ct_eq(&self, other: &Self) -> subtle::Choice {
+        CtEq::ct_eq(self, other).into()
+    }
+}
+
+#[cfg(feature = "subtle")]
+impl<const LIMBS: usize> subtle::ConditionallySelectable for PrimeParams<LIMBS> {
+    fn conditional_assign(&mut self, src: &Self, choice: subtle::Choice) {
+        self.ct_assign(src, choice.into());
+    }
+
+    fn conditional_select(a: &Self, b: &Self, choice: subtle::Choice) -> Self {
+        a.ct_select(b, choice.into())
+    }
+}
+
+#[allow(clippy::unwrap_in_result)]
+const fn find_primitive_root<const LIMBS: usize>(
+    p: &Odd<Uint<LIMBS>>,
+) -> Option<(NonZeroU32, Uint<LIMBS>)> {
+    // A primitive root exists iff p is 1, 2, 4, q^k or 2q^k, k > 0, q is an odd prime.
+    // Find a quadratic non-residue (primitive roots are non-residue for powers of a prime)
+    let mut g = NonZeroU32::new(2u32).expect("ensured non-zero");
+    loop {
+        // Either the modulus is prime and g is quadratic non-residue, or
+        // the modulus is composite.
+        let g_uint = Uint::from_u32(g.get());
+        match g_uint.jacobi_symbol_vartime(p) as i8 {
+            -1 => {
+                break Some((g, g_uint));
+            }
+            0 => {
+                // Modulus is composite
+                return None;
+            }
+            _ => {
+                let Some(g2) = g.checked_add(1) else {
+                    return None;
+                };
+                g = g2;
+            }
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::PrimeParams;
+    use crate::{Choice, CtEq, CtSelect, Odd, U128, modular::MontyParams};
+
+    #[test]
+    fn check_expected() {
+        let monty_params =
+            MontyParams::new_vartime(Odd::<U128>::from_be_hex("e38af050d74b8567f73c8713cbc7bc47"));
+        let prime_params = PrimeParams::new_vartime(&monty_params).expect("failed creating params");
+        assert_eq!(prime_params.s.get(), 1);
+        assert_eq!(prime_params.generator.get(), 5);
+    }
+
+    #[test]
+    fn check_non_prime() {
+        let monty_params =
+            MontyParams::new_vartime(Odd::<U128>::from_be_hex("e38af050d74b8567f73c8713cbc7bc01"));
+        assert!(PrimeParams::new_vartime(&monty_params).is_none());
+    }
+
+    #[test]
+    fn check_equality() {
+        let monty_params_1 =
+            MontyParams::new_vartime(Odd::<U128>::from_be_hex("e38af050d74b8567f73c8713cbc7bc47"));
+        let prime_params_1 =
+            PrimeParams::new_vartime(&monty_params_1).expect("failed creating params");
+
+        let monty_params_2 =
+            MontyParams::new_vartime(Odd::<U128>::from_be_hex("f2799d643ab7ff983437c3a86cdb1beb"));
+        let prime_params_2 =
+            PrimeParams::new_vartime(&monty_params_2).expect("failed creating params");
+
+        assert!(CtEq::ct_eq(&prime_params_1, &prime_params_1).to_bool_vartime());
+        #[cfg(feature = "subtle")]
+        assert!(bool::from(subtle::ConstantTimeEq::ct_eq(
+            &prime_params_1,
+            &prime_params_1
+        )));
+
+        assert!(CtEq::ct_ne(&prime_params_1, &prime_params_2).to_bool_vartime());
+        #[cfg(feature = "subtle")]
+        assert!(bool::from(subtle::ConstantTimeEq::ct_ne(
+            &prime_params_1,
+            &prime_params_2
+        )));
+
+        assert!(
+            CtEq::ct_eq(
+                &CtSelect::ct_select(&prime_params_1, &prime_params_2, Choice::FALSE),
+                &prime_params_1,
+            )
+            .to_bool_vartime()
+        );
+        #[cfg(feature = "subtle")]
+        assert!(
+            CtEq::ct_eq(
+                &subtle::ConditionallySelectable::conditional_select(
+                    &prime_params_1,
+                    &prime_params_2,
+                    subtle::Choice::from(0u8)
+                ),
+                &prime_params_1,
+            )
+            .to_bool_vartime()
+        );
+
+        assert!(
+            CtEq::ct_eq(
+                &CtSelect::ct_select(&prime_params_1, &prime_params_2, Choice::TRUE),
+                &prime_params_2,
+            )
+            .to_bool_vartime()
+        );
+        #[cfg(feature = "subtle")]
+        assert!(
+            CtEq::ct_eq(
+                &subtle::ConditionallySelectable::conditional_select(
+                    &prime_params_1,
+                    &prime_params_2,
+                    subtle::Choice::from(1u8)
+                ),
+                &prime_params_2,
+            )
+            .to_bool_vartime()
+        );
+    }
+}
diff --git a/src/modular/sqrt.rs b/src/modular/sqrt.rs
new file mode 100644
index 00000000..3e40eaa3
--- /dev/null
+++ b/src/modular/sqrt.rs
@@ -0,0 +1,295 @@
+use crate::{
+    Choice, CtOption, Uint,
+    modular::{FixedMontyForm, FixedMontyParams, prime_params::PrimeParams},
+};
+
+/// Compute a modular square root (if it exists) given [`MontyParams`]
+/// and [`PrimeParams`] corresponding to `monty_value`, in Montgomery form.
+#[must_use]
+pub const fn sqrt_montgomery_form<const LIMBS: usize>(
+    monty_value: &Uint<LIMBS>,
+    monty_params: &FixedMontyParams<LIMBS>,
+    prime_params: &PrimeParams<LIMBS>,
+) -> CtOption<Uint<LIMBS>> {
+    let value = FixedMontyForm::from_montgomery(*monty_value, monty_params);
+    let b = value.pow_vartime(&prime_params.sqrt_exp);
+
+    // Constant-time versions of modular square root algorithms based on:
+    // Koo, N., Cho, G.H. and Kwon, S. (2013), "Square root algorithm in 𝔽q for q ≡ 2s + 1 (mod 2s+1)".
+    // Electron. Lett., 49: 467-469. https://doi.org/10.1049/el.2012.4239
+
+    let x = match prime_params.s.get() {
+        1 => {
+            // Shanks algorithm: sqrt = x^((p+1)/4) = x^(t+1)
+            b
+        }
+        2 => {
+            // Algorithm 3: p = 5 mod 8 (Atkins variant)
+            let ru = FixedMontyForm::from_montgomery(prime_params.monty_root_unity, monty_params);
+            let cb = value.mul(&b);
+            let zeta = cb.mul(&b);
+            let is_one = Uint::eq(zeta.as_montgomery(), monty_params.one());
+            monty_select(&cb.mul(&ru), &cb, is_one)
+        }
+        3 => {
+            // Algorithm 4: p = 9 mod 16
+            let ru = FixedMontyForm::from_montgomery(prime_params.monty_root_unity, monty_params);
+            let ru_2 =
+                FixedMontyForm::from_montgomery(prime_params.monty_root_unity_p2, monty_params);
+            let ru_3 = ru.mul(&ru_2);
+            let cb = value.mul(&b);
+            let zeta = cb.mul(&b);
+
+            let mut m = monty_select(
+                &ru,
+                &FixedMontyForm::one(ru.params()),
+                Uint::eq(zeta.as_montgomery(), monty_params.one()),
+            );
+            // m = ru^2 if zeta = -1
+            m = monty_select(
+                &m,
+                &ru_2,
+                Uint::eq(zeta.neg().as_montgomery(), monty_params.one()),
+            );
+            // m = ru^3 if zeta = ru^2
+            m = monty_select(&m, &ru_3, monty_eq(&zeta, &ru_2));
+
+            cb.mul(&m)
+        }
+        4 => {
+            // Algorithm 5: p = 17 mod 32
+            let ru = FixedMontyForm::from_montgomery(prime_params.monty_root_unity, monty_params);
+            let ru_2 =
+                FixedMontyForm::from_montgomery(prime_params.monty_root_unity_p2, monty_params);
+            let ru_4 = ru_2.square();
+            let ru_6 = ru_2.mul(&ru_4);
+            let cb = value.mul(&b);
+            let zeta = cb.mul(&b);
+
+            let neg_zeta = zeta.neg();
+            let zeta_b = monty_eq(&zeta, &ru_2);
+            let neg_zeta_b = monty_eq(&neg_zeta, &ru_2);
+            let zeta_d = monty_eq(&zeta, &ru_6);
+
+            // m = B if -zeta in (B, C), else 1
+            let mut m = monty_select(
+                &FixedMontyForm::one(ru.params()),
+                &ru_2,
+                neg_zeta_b.or(monty_eq(&neg_zeta, &ru_4)),
+            );
+            // m = C if zeta in (-1, D)
+            m = monty_select(
+                &m,
+                &ru_4,
+                Uint::eq(neg_zeta.as_montgomery(), monty_params.one()).or(zeta_d),
+            );
+            // m = D if zeta in (B, C)
+            m = monty_select(&m, &ru_6, zeta_b.or(monty_eq(&zeta, &ru_4)));
+            // m = m•ru if zeta or -zeta in (B, D)
+            m = monty_select(
+                &m,
+                &m.mul(&ru),
+                zeta_b
+                    .or(zeta_d)
+                    .or(neg_zeta_b)
+                    .or(monty_eq(&neg_zeta, &ru_6)),
+            );
+
+            cb.mul(&m)
+        }
+        _ => {
+            // Tonelli-Shanks
+            let mut x = value.mul(&b);
+            let mut d = x.mul(&b);
+            let mut z =
+                FixedMontyForm::from_montgomery(prime_params.monty_root_unity, monty_params);
+            let mut v = prime_params.s.get();
+            let mut max_v = v;
+
+            while max_v >= 1 {
+                let mut k = 1;
+                let mut tmp = d.square();
+                let mut j_less_than_v = Choice::TRUE;
+
+                let mut j = 2;
+                while j < max_v {
+                    let tmp_is_one = Uint::eq(tmp.as_montgomery(), monty_params.one());
+                    let squared = monty_select(&tmp, &z, tmp_is_one).square();
+                    tmp = monty_select(&squared, &tmp, tmp_is_one);
+                    j_less_than_v = j_less_than_v.and(Choice::from_u32_eq(j, v).not());
+                    z = monty_select(&z, &squared, tmp_is_one.and(j_less_than_v));
+                    k = tmp_is_one.select_u32(j, k);
+                    j += 1;
+                }
+
+                let b_is_one = Uint::eq(d.as_montgomery(), monty_params.one());
+                x = monty_select(&x.mul(&z), &x, b_is_one);
+                z = z.square();
+                d = d.mul(&z);
+                v = k;
+                max_v -= 1;
+            }
+
+            x
+        }
+    };
+
+    CtOption::new(x.to_montgomery(), monty_eq(&x.square(), &value))
+}
+
+const fn monty_eq<const LIMBS: usize>(
+    a: &FixedMontyForm<LIMBS>,
+    b: &FixedMontyForm<LIMBS>,
+) -> Choice {
+    Uint::eq(a.as_montgomery(), b.as_montgomery())
+}
+
+const fn monty_select<const LIMBS: usize>(
+    a: &FixedMontyForm<LIMBS>,
+    b: &FixedMontyForm<LIMBS>,
+    c: Choice,
+) -> FixedMontyForm<LIMBS> {
+    FixedMontyForm::from_montgomery(
+        Uint::select(a.as_montgomery(), b.as_montgomery(), c),
+        a.params(),
+    )
+}
+
+#[cfg(test)]
+mod tests {
+    use super::sqrt_montgomery_form;
+    use crate::{
+        Odd, U256, U576, Uint,
+        modular::{FixedMontyForm, FixedMontyParams, PrimeParams},
+    };
+
+    fn root_of_unity<const LIMBS: usize>(
+        monty_params: &FixedMontyParams<LIMBS>,
+        prime_params: &PrimeParams<LIMBS>,
+    ) -> Uint<LIMBS> {
+        FixedMontyForm::from_montgomery(prime_params.monty_root_unity, monty_params).retrieve()
+    }
+
+    fn test_monty_sqrt<const LIMBS: usize>(
+        monty_params: FixedMontyParams<LIMBS>,
+        prime_params: PrimeParams<LIMBS>,
+    ) {
+        let modulus = monty_params.modulus.get();
+        let rounds = if cfg!(miri) { 1..=2 } else { 0..=256 };
+        for i in rounds {
+            let s = i * i;
+            let s_monty = FixedMontyForm::new(&Uint::from_u32(s), &monty_params);
+            let rt_monty =
+                sqrt_montgomery_form(s_monty.as_montgomery(), &monty_params, &prime_params)
+                    .expect("no sqrt found");
+            let rt = FixedMontyForm::from_montgomery(rt_monty, &monty_params).retrieve();
+            let i = Uint::from_u32(i);
+            assert!(
+                Uint::eq(&rt, &i)
+                    .or(Uint::eq(&rt, &modulus.wrapping_sub(&i)))
+                    .to_bool_vartime()
+            );
+        }
+
+        // generator must be non-residue
+        let generator = Uint::from_u32(prime_params.generator.get());
+        let gen_monty = FixedMontyForm::new(&generator, &monty_params);
+        assert!(
+            sqrt_montgomery_form(gen_monty.as_montgomery(), &monty_params, &prime_params)
+                .is_none()
+                .to_bool_vartime()
+        );
+    }
+
+    #[test]
+    fn mod_sqrt_s_1() {
+        // p = 3 mod 4, s = 1
+        // P-256 field modulus
+        let monty_params = FixedMontyParams::new_vartime(Odd::<U256>::from_be_hex(
+            "ffffffff00000001000000000000000000000000ffffffffffffffffffffffff",
+        ));
+        let prime_params = PrimeParams::new_vartime(&monty_params).expect("failed creating params");
+        assert_eq!(prime_params.s.get(), 1);
+        assert_eq!(prime_params.generator.get(), 3);
+        assert_eq!(
+            root_of_unity(&monty_params, &prime_params),
+            U256::from_be_hex("FFFFFFFF00000001000000000000000000000000FFFFFFFFFFFFFFFFFFFFFFFE")
+        );
+
+        test_monty_sqrt(monty_params, prime_params);
+    }
+
+    #[test]
+    fn mod_sqrt_s_2() {
+        // p = 5 mod 8, s = 2
+        // ed25519 base field
+        let monty_params = FixedMontyParams::new_vartime(Odd::<U256>::from_be_hex(
+            "7fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffed",
+        ));
+        let prime_params = PrimeParams::new_vartime(&monty_params).expect("failed creating params");
+        assert_eq!(prime_params.s.get(), 2);
+        assert_eq!(prime_params.generator.get(), 2);
+        assert_eq!(
+            root_of_unity(&monty_params, &prime_params),
+            U256::from_be_hex("2B8324804FC1DF0B2B4D00993DFBD7A72F431806AD2FE478C4EE1B274A0EA0B0")
+        );
+
+        test_monty_sqrt(monty_params, prime_params);
+    }
+
+    #[test]
+    fn mod_sqrt_s_3() {
+        // p = 9 mod 16, s = 3
+        // brainpoolP384 scalar field
+        let monty_params = FixedMontyParams::new_vartime(Odd::<U576>::from_be_hex(
+            "00000000000001fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffa51868783bf2f966b7fcc0148f709a5d03bb5c9b8899c47aebb6fb71e91386409",
+        ));
+        let prime_params = PrimeParams::new_vartime(&monty_params).expect("failed creating params");
+        assert_eq!(prime_params.s.get(), 3);
+        assert_eq!(prime_params.generator.get(), 3);
+        assert_eq!(
+            root_of_unity(&monty_params, &prime_params),
+            U576::from_be_hex(
+                "000000000000009a0a650d44b28c17f3d708ad2fa8c4fbc7e6000d7c12dafa92fcc5673a3055276d535f79ff391dcdbcd998b7836647d3a72472b3da861ac810a7f9c7b7b63e2205"
+            )
+        );
+
+        test_monty_sqrt(monty_params, prime_params);
+    }
+
+    #[test]
+    fn mod_sqrt_s_4() {
+        // p = 17 mod 32, s = 4
+        // P-256 scalar field
+        let monty_params = FixedMontyParams::new_vartime(Odd::<U256>::from_be_hex(
+            "ffffffff00000000ffffffffffffffffbce6faada7179e84f3b9cac2fc632551",
+        ));
+        let prime_params = PrimeParams::new_vartime(&monty_params).expect("failed creating params");
+        assert_eq!(prime_params.s.get(), 4);
+        assert_eq!(prime_params.generator.get(), 7);
+        assert_eq!(
+            root_of_unity(&monty_params, &prime_params),
+            U256::from_be_hex("ffc97f062a770992ba807ace842a3dfc1546cad004378daf0592d7fbb41e6602")
+        );
+
+        test_monty_sqrt(monty_params, prime_params);
+    }
+
+    #[test]
+    fn mod_sqrt_s_6() {
+        // s = 6
+        // K-256 scalar field
+        let monty_params = FixedMontyParams::new_vartime(Odd::<U256>::from_be_hex(
+            "FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEBAAEDCE6AF48A03BBFD25E8CD0364141",
+        ));
+        let prime_params = PrimeParams::new_vartime(&monty_params).expect("failed creating params");
+        assert_eq!(prime_params.s.get(), 6);
+        assert_eq!(prime_params.generator.get(), 5);
+        assert_eq!(
+            root_of_unity(&monty_params, &prime_params),
+            U256::from_be_hex("0D1F8EAB98DCD1ACA7DC810E065710CBB96E9ABEBBE451FA15B4F83D2D2AD232")
+        );
+
+        test_monty_sqrt(monty_params, prime_params);
+    }
+}