diff --git a/Cargo.lock b/Cargo.lock
index f6fd77a28..8d9ca6a3f 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -6789,10 +6789,13 @@ dependencies = [
 name = "openvm-accelerators"
 version = "0.4.0"
 dependencies = [
+ "aurora-engine-modexp",
  "hex-literal",
  "k256 0.13.4 (registry+https://github.com/rust-lang/crates.io-index)",
  "k256 0.13.4 (git+https://github.com/openvm-org/openvm.git?branch=develop-v2.1.0)",
+ "openvm-ecc-guest",
  "openvm-keccak256",
+ "openvm-pairing",
  "openvm-sha2",
  "p256 0.13.2 (git+https://github.com/openvm-org/openvm.git?branch=develop-v2.1.0)",
  "ripemd",
diff --git a/crates/accelerators/Cargo.toml b/crates/accelerators/Cargo.toml
index 19648018d..e89886d92 100644
--- a/crates/accelerators/Cargo.toml
+++ b/crates/accelerators/Cargo.toml
@@ -11,11 +11,14 @@ workspace = true
 
 [dependencies]
 # openvm
+openvm-ecc-guest.workspace = true
 openvm-p256.workspace = true
+openvm-pairing = { workspace = true, features = ["bn254"] }
 openvm-keccak256.workspace = true
 openvm-sha2.workspace = true
 
 # crypto
+aurora-engine-modexp = { version = "1.2.0", default-features = false }
 ripemd = { version = "0.1.3", default-features = false }
 
 # The OpenVM-accelerated k256 fork; its ECDSA recovery relies on zkVM hints
diff --git a/crates/accelerators/src/ffi/mod.rs b/crates/accelerators/src/ffi/mod.rs
index 3f82ef4ea..e351d6812 100644
--- a/crates/accelerators/src/ffi/mod.rs
+++ b/crates/accelerators/src/ffi/mod.rs
@@ -7,7 +7,9 @@
 mod blake2;
 mod ecdsa;
 mod hash;
+mod modexp;
 
 pub use blake2::*;
 pub use ecdsa::*;
 pub use hash::*;
+pub use modexp::*;
diff --git a/crates/accelerators/src/ffi/modexp.rs b/crates/accelerators/src/ffi/modexp.rs
new file mode 100644
index 000000000..3733d987f
--- /dev/null
+++ b/crates/accelerators/src/ffi/modexp.rs
@@ -0,0 +1,48 @@
+//! C ABI for modular exponentiation.
+
+use crate::{ops, types::ZkvmStatus};
+
+/// Compute `base[..base_len] ^ exp[..exp_len] % modulus[..mod_len]` into
+/// `output`, which receives exactly `mod_len` bytes, left-padded with zeros.
+///
+/// Returns [`ZkvmStatus::Fail`] if any pointer is NULL while its length is
+/// non-zero; a NULL pointer with a zero length is the empty input.
+///
+/// # Safety
+///
+/// - `base`, `exp` and `modulus`, if non-NULL, must be valid for reads of their respective lengths.
+/// - `output`, if non-NULL, must be valid for writes of `mod_len` bytes.
+#[unsafe(no_mangle)]
+pub unsafe extern "C" fn zkvm_modexp(
+    base: *const u8,
+    base_len: usize,
+    exp: *const u8,
+    exp_len: usize,
+    modulus: *const u8,
+    mod_len: usize,
+    output: *mut u8,
+) -> ZkvmStatus {
+    if (base.is_null() && base_len != 0) ||
+        (exp.is_null() && exp_len != 0) ||
+        (modulus.is_null() && mod_len != 0) ||
+        (output.is_null() && mod_len != 0)
+    {
+        return ZkvmStatus::Fail;
+    }
+    // SAFETY: non-NULL checked above; validity is guaranteed by the caller.
+    let base =
+        if base_len == 0 { &[] } else { unsafe { core::slice::from_raw_parts(base, base_len) } };
+    // SAFETY: non-NULL checked above; validity is guaranteed by the caller.
+    let exp = if exp_len == 0 { &[] } else { unsafe { core::slice::from_raw_parts(exp, exp_len) } };
+    // SAFETY: non-NULL checked above; validity is guaranteed by the caller.
+    let modulus =
+        if mod_len == 0 { &[] } else { unsafe { core::slice::from_raw_parts(modulus, mod_len) } };
+    // SAFETY: non-NULL checked above; validity is guaranteed by the caller.
+    let output = if mod_len == 0 {
+        &mut [][..]
+    } else {
+        unsafe { core::slice::from_raw_parts_mut(output, mod_len) }
+    };
+    ops::modexp(base, exp, modulus, output);
+    ZkvmStatus::Ok
+}
diff --git a/crates/accelerators/src/lib.rs b/crates/accelerators/src/lib.rs
index 10a85ec01..e851bff8e 100644
--- a/crates/accelerators/src/lib.rs
+++ b/crates/accelerators/src/lib.rs
@@ -2,6 +2,8 @@
 
 #![cfg_attr(not(feature = "std"), no_std)]
 
+extern crate alloc;
+
 #[cfg(feature = "ffi")]
 pub mod ffi;
 pub mod ops;
diff --git a/crates/accelerators/src/ops/mod.rs b/crates/accelerators/src/ops/mod.rs
index f269b7319..4813df2b2 100644
--- a/crates/accelerators/src/ops/mod.rs
+++ b/crates/accelerators/src/ops/mod.rs
@@ -7,10 +7,12 @@
 mod blake2;
 mod ecdsa;
 mod hash;
+mod modexp;
 
 pub use blake2::blake2f;
 pub use ecdsa::{secp256k1_ecrecover, secp256k1_verify, secp256r1_verify};
 pub use hash::{keccak256, ripemd160, sha256};
+pub use modexp::modexp;
 
 #[derive(Clone, Copy, Debug, PartialEq, Eq)]
 pub enum Error {
diff --git a/crates/accelerators/src/ops/modexp.rs b/crates/accelerators/src/ops/modexp.rs
new file mode 100644
index 000000000..5652d7362
--- /dev/null
+++ b/crates/accelerators/src/ops/modexp.rs
@@ -0,0 +1,139 @@
+//! Modular exponentiation with a BN254-Fr fast path.
+
+use alloc::{vec, vec::Vec};
+
+use openvm_ecc_guest::algebra::{ExpBytes, IntMod, Reduce};
+use openvm_pairing::bn254 as bn;
+
+/// The number of bytes needed to represent an element of BN254's scalar
+/// field Fr.
+const BN_SCALAR_LEN: usize = 32;
+
+/// Compute `base^exp % modulus` into `output`, left-padded with zeros.
+///
+/// # Panics
+///
+/// Panics if `output.len() != modulus.len()`.
+pub fn modexp(base: &[u8], exp: &[u8], modulus: &[u8], output: &mut [u8]) {
+    assert_eq!(output.len(), modulus.len(), "output must be exactly modulus-sized");
+
+    let result = if is_bn254_fr(modulus) {
+        accelerated_modexp_bn254_fr(base, exp)
+    } else {
+        aurora_engine_modexp::modexp(base, exp, modulus)
+    };
+
+    // The result is numerically reduced, but its byte representation may be
+    // shorter or longer (leading zeros); right-align it.
+    if result.len() >= output.len() {
+        output.copy_from_slice(&result[result.len() - output.len()..]);
+    } else {
+        let pad = output.len() - result.len();
+        output[..pad].fill(0);
+        output[pad..].copy_from_slice(&result);
+    }
+}
+
+/// Returns true if the modulus (big-endian, possibly with leading zeros) equals BN254 Fr.
+fn is_bn254_fr(modulus: &[u8]) -> bool {
+    // Strip leading zeros
+    let stripped = match modulus.iter().position(|&b| b != 0) {
+        Some(i) => &modulus[i..],
+        None => return false, // all zeros
+    };
+    // bn::Scalar::MODULUS is little-endian; compare against reversed input
+    stripped.len() == BN_SCALAR_LEN && stripped.iter().rev().eq(bn::Scalar::MODULUS.as_ref().iter())
+}
+
+/// Accelerated modexp for BN254 Fr using field arithmetic intrinsics.
+fn accelerated_modexp_bn254_fr(base: &[u8], exp: &[u8]) -> Vec<u8> {
+    // OpenVM's field reduction requires inputs to be aligned to the field byte size.
+    let padded_len = base.len().next_multiple_of(BN_SCALAR_LEN).max(BN_SCALAR_LEN);
+    let mut padded = vec![0u8; padded_len];
+    padded[padded_len - base.len()..].copy_from_slice(base);
+    let base_fr = bn::Scalar::reduce_be_bytes(&padded);
+
+    base_fr.exp_bytes(true, exp).to_be_bytes().as_ref().to_vec()
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    /// BN254 Fr modulus in big-endian bytes
+    fn bn254_fr_modulus_be() -> Vec<u8> {
+        bn::Scalar::MODULUS.as_ref().iter().rev().copied().collect()
+    }
+
+    /// Helper: run the accelerated path and compare against the aurora
+    /// reference. The accelerated path always returns BN_SCALAR_LEN bytes,
+    /// so the reference output is left-padded to match.
+    fn check(base: &[u8], exp: &[u8]) {
+        let modulus = bn254_fr_modulus_be();
+        let expected = aurora_engine_modexp::modexp(base, exp, &modulus);
+        let actual = accelerated_modexp_bn254_fr(base, exp);
+        let mut expected_padded = vec![0u8; BN_SCALAR_LEN];
+        let offset = BN_SCALAR_LEN - expected.len();
+        expected_padded[offset..].copy_from_slice(&expected);
+        assert_eq!(actual, expected_padded, "base={base:?}, exp={exp:?}");
+    }
+
+    #[test]
+    fn test_is_bn254_fr() {
+        // Exact modulus
+        assert!(is_bn254_fr(&bn254_fr_modulus_be()));
+
+        // With leading zeros
+        let mut padded = vec![0u8; 10];
+        padded.extend_from_slice(&bn254_fr_modulus_be());
+        assert!(is_bn254_fr(&padded));
+
+        // All zeros → false
+        assert!(!is_bn254_fr(&[0u8; 32]));
+
+        // Wrong modulus (flip last bit)
+        let mut m = bn254_fr_modulus_be();
+        *m.last_mut().unwrap() ^= 1;
+        assert!(!is_bn254_fr(&m));
+    }
+
+    #[test]
+    fn test_accelerated_modexp_bn254_fr() {
+        // --- short base (<=32 bytes), value < modulus ---
+        check(&[3], &[5]); // 3^5 mod Fr
+        check(&[0], &[5]); // 0^5 = 0
+        check(&[3], &[0]); // 3^0 = 1
+        check(&[0], &[0]); // 0^0 = 1 by convention
+        check(&[], &[]); // empty inputs
+        check(&[0, 0, 0, 3], &[5]); // leading zeros in base
+
+        // --- short base, value >= modulus (triggers the reduce fallback) ---
+        let m = bn254_fr_modulus_be();
+        check(&m, &[1]); // Fr mod Fr = 0, so 0^1 = 0
+        let mut m_plus_1 = m.clone();
+        *m_plus_1.last_mut().unwrap() += 1;
+        check(&m_plus_1, &[2]); // (Fr+1)^2 mod Fr = 1
+        check(&[0xff; 32], &[1]); // max 256-bit value, >= modulus
+
+        // --- large base (> 32 bytes, reduce_be_bytes path) ---
+        check(&[0xab; 64], &[3]); // aligned (multiple of 32)
+        check(&[0x42; 100], &[2]); // unaligned (tests the padding)
+        check(&[0xab; 64], &[0xff; 32]); // large base + large exponent
+
+        // --- larger exponents ---
+        check(&[2], &[0xff; 32]); // 2^(2^256-1) mod Fr
+        check(&[2], &[0, 0, 0, 5]); // leading zeros in exponent
+        check(&[3], &[0xab; 64]); // exponent > 32 bytes
+
+        // --- same value through both base-parsing code paths ---
+        let base_32 = [0xab; 32];
+        let mut base_33 = vec![0u8];
+        base_33.extend_from_slice(&base_32);
+        let exp = &[7];
+        assert_eq!(
+            accelerated_modexp_bn254_fr(&base_32, exp),
+            accelerated_modexp_bn254_fr(&base_33, exp),
+            "33-byte base with leading zero must match 32-byte base"
+        );
+    }
+}
diff --git a/crates/accelerators/tests/conformance/main.rs b/crates/accelerators/tests/conformance/main.rs
index d8bb24d1e..af50b7471 100644
--- a/crates/accelerators/tests/conformance/main.rs
+++ b/crates/accelerators/tests/conformance/main.rs
@@ -6,3 +6,4 @@
 mod blake2;
 mod ecdsa;
 mod hash;
+mod modexp;
diff --git a/crates/accelerators/tests/conformance/modexp.rs b/crates/accelerators/tests/conformance/modexp.rs
new file mode 100644
index 000000000..760015463
--- /dev/null
+++ b/crates/accelerators/tests/conformance/modexp.rs
@@ -0,0 +1,107 @@
+//! Modexp conformance vectors.
+
+use hex_literal::hex;
+use openvm_accelerators::{ffi::zkvm_modexp, ops::modexp, types::ZkvmStatus};
+
+/// BN254 Fr (the scalar field) modulus, big-endian. Not to be confused with
+/// the base field prime, which shares the leading bytes.
+const BN254_FR: [u8; 32] = hex!("30644e72e131a029b85045b68181585d2833e84879b9709143e1f593f0000001");
+
+#[test]
+fn modexp_small() {
+    // 3^5 mod 7 = 5
+    let mut output = [0xffu8; 1];
+    modexp(&[3], &[5], &[7], &mut output);
+    assert_eq!(output, [5]);
+
+    // Output is left-padded to the modulus length.
+    let mut output = [0xffu8; 2];
+    modexp(&[3], &[5], &[0, 7], &mut output);
+    assert_eq!(output, [0, 5]);
+
+    // A zero-length modulus writes nothing.
+    modexp(&[3], &[5], &[], &mut []);
+}
+
+#[test]
+fn modexp_matches_reference() {
+    // The BN254-Fr accelerated path, compared right-aligned against the
+    // aurora reference.
+    let mut output = [0u8; 32];
+    modexp(&[0xab; 32], &[0x07], &BN254_FR, &mut output);
+    let reference = aurora_engine_modexp::modexp(&[0xab; 32], &[0x07], &BN254_FR);
+    assert_eq!(output[32 - reference.len()..], reference[..]);
+
+    // The generic path with a non-special modulus.
+    let modulus = [0xef; 24];
+    let mut output = [0u8; 24];
+    modexp(&[0x12; 40], &[0x34; 3], &modulus, &mut output);
+    let reference = aurora_engine_modexp::modexp(&[0x12; 40], &[0x34; 3], &modulus);
+    assert_eq!(output[24 - reference.len()..], reference[..]);
+}
+
+#[test]
+fn zkvm_modexp_smoke() {
+    // 3^5 mod 7 = 5
+    let base = [3u8];
+    let exp = [5u8];
+    let modulus = [7u8];
+    let mut output = [0xffu8; 1];
+    let status = unsafe {
+        zkvm_modexp(base.as_ptr(), 1, exp.as_ptr(), 1, modulus.as_ptr(), 1, output.as_mut_ptr())
+    };
+    assert_eq!(status, ZkvmStatus::Ok);
+    assert_eq!(output, [5]);
+}
+
+#[test]
+fn zkvm_modexp_null_pointers() {
+    let base = [3u8];
+    let exp = [5u8];
+    let modulus = [7u8];
+    let mut output = [0xffu8; 1];
+
+    // NULL base and exp with zero lengths are empty inputs: 0^0 mod 7 = 1.
+    let status = unsafe {
+        zkvm_modexp(
+            core::ptr::null(),
+            0,
+            core::ptr::null(),
+            0,
+            modulus.as_ptr(),
+            1,
+            output.as_mut_ptr(),
+        )
+    };
+    assert_eq!(status, ZkvmStatus::Ok);
+    assert_eq!(output, [1]);
+
+    // A NULL pointer with a non-zero length fails.
+    let status = unsafe {
+        zkvm_modexp(core::ptr::null(), 1, exp.as_ptr(), 1, modulus.as_ptr(), 1, output.as_mut_ptr())
+    };
+    assert_eq!(status, ZkvmStatus::Fail);
+
+    let status = unsafe {
+        zkvm_modexp(
+            base.as_ptr(),
+            1,
+            core::ptr::null(),
+            1,
+            modulus.as_ptr(),
+            1,
+            output.as_mut_ptr(),
+        )
+    };
+    assert_eq!(status, ZkvmStatus::Fail);
+
+    let status = unsafe {
+        zkvm_modexp(base.as_ptr(), 1, exp.as_ptr(), 1, core::ptr::null(), 1, output.as_mut_ptr())
+    };
+    assert_eq!(status, ZkvmStatus::Fail);
+
+    let status = unsafe {
+        zkvm_modexp(base.as_ptr(), 1, exp.as_ptr(), 1, modulus.as_ptr(), 1, core::ptr::null_mut())
+    };
+    assert_eq!(status, ZkvmStatus::Fail);
+}