diff --git a/Cargo.lock b/Cargo.lock
index 7559226..58553de 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -58,6 +58,7 @@ version = "0.1.0"
 dependencies = [
  "curve",
  "field",
+ "rand",
 ]
 
 [[package]]
diff --git a/isogeny/Cargo.toml b/isogeny/Cargo.toml
index 12e5bc1..c5bc4fe 100644
--- a/isogeny/Cargo.toml
+++ b/isogeny/Cargo.toml
@@ -6,3 +6,6 @@ edition = "2024"
 [dependencies]
 field = { path = "../field" }
 curve = { path = "../curve" }
+
+[dev-dependencies]
+rand = "0.10.1"
diff --git a/isogeny/examples/csidh.rs b/isogeny/examples/csidh.rs
new file mode 100644
index 0000000..5cfb3c2
--- /dev/null
+++ b/isogeny/examples/csidh.rs
@@ -0,0 +1,43 @@
+use curve::Curve;
+use isogeny::csidh;
+use rand::{RngExt, rng};
+
+const KEY_BOUND: usize = 74;
+
+fn random_private_key<R: RngExt>(rng: &mut R) -> [i8; csidh::PRIMES.len()] {
+    let mut key = [0i8; csidh::PRIMES.len()];
+    for slot in key.iter_mut().take(KEY_BOUND) {
+        *slot = (rng.random::<u8>() % 3) as i8 - 1;
+    }
+    key
+}
+
+fn main() {
+    let mut rng = rng();
+    let a_priv = random_private_key(&mut rng);
+    let b_priv = random_private_key(&mut rng);
+    let ells: Vec<u64> = csidh::PRIMES.iter().take(KEY_BOUND).copied().collect();
+    println!("Small primes in play: {:?}", ells);
+    println!("Alice's private:      {:?}", &a_priv[..KEY_BOUND]);
+    println!("Bob's private:        {:?}", &b_priv[..KEY_BOUND]);
+
+    let start = csidh::base_curve();
+
+    println!("\nDeriving public curves...");
+    let t = std::time::Instant::now();
+    let a_pub = csidh::action(&start, &a_priv, &mut || rng.random::<u64>());
+    let b_pub = csidh::action(&start, &b_priv, &mut || rng.random::<u64>());
+    println!("done in {:?}", t.elapsed());
+
+    println!("\nDeriving shared secrets...");
+    let t = std::time::Instant::now();
+    let a_shared = csidh::action(&b_pub, &a_priv, &mut || rng.random::<u64>());
+    let b_shared = csidh::action(&a_pub, &b_priv, &mut || rng.random::<u64>());
+    println!("done in {:?}", t.elapsed());
+
+    let j_a = a_shared.j_invariant();
+    let j_b = b_shared.j_invariant();
+    assert_eq!(j_a, j_b, "shared j-invariants disagreed");
+    println!("\nShared j-invariant: {:?}", j_a);
+    println!("CSIDH shared secret established.");
+}
diff --git a/isogeny/src/lib.rs b/isogeny/src/lib.rs
index 8b13789..d5dafe7 100644
--- a/isogeny/src/lib.rs
+++ b/isogeny/src/lib.rs
@@ -1 +1,366 @@
+use curve::{Curve, Point, PointRepresentation};
+use field::{Field, FieldElement};
 
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub struct ShortWeierstrass<F: Field> {
+    pub a: FieldElement<F>,
+    pub b: FieldElement<F>,
+}
+
+impl<F: Field> Curve<F> for ShortWeierstrass<F> {
+    fn a(&self) -> FieldElement<F> {
+        self.a
+    }
+
+    fn b(&self) -> FieldElement<F> {
+        self.b
+    }
+}
+
+#[derive(Debug, Clone)]
+pub struct Isogeny<F: Field> {
+    domain: ShortWeierstrass<F>,
+    codomain: ShortWeierstrass<F>,
+    s: Vec<(FieldElement<F>, FieldElement<F>, FieldElement<F>)>,
+}
+
+impl<F: Field> Isogeny<F> {
+    pub fn velu<C: Curve<F>>(domain: &C, kernel: Vec<Point<F>>) -> Self {
+        let a = domain.a();
+        let b = domain.b();
+
+        let mut reps: Vec<(FieldElement<F>, FieldElement<F>)> = Vec::new();
+        for p in kernel {
+            if let PointRepresentation::Affine { x, y } = p.representation() {
+                if !reps.iter().any(|(rx, _)| *rx == x) {
+                    reps.push((x, y));
+                }
+            }
+        }
+
+        let mut v_sum = FieldElement::zero();
+        let mut w_sum = FieldElement::zero();
+        let mut s = Vec::with_capacity(reps.len());
+        for (x_q, y_q) in reps {
+            let g_x = FieldElement::three() * x_q * x_q + a;
+            let v_q = if y_q == FieldElement::zero() {
+                g_x
+            } else {
+                FieldElement::two() * g_x
+            };
+            let u_q = FieldElement::from_u64(4) * y_q * y_q;
+            v_sum += v_q;
+            w_sum += u_q + x_q * v_q;
+            s.push((x_q, v_q, u_q));
+        }
+
+        let codomain = ShortWeierstrass {
+            a: a - FieldElement::from_u64(5) * v_sum,
+            b: b - FieldElement::from_u64(7) * w_sum,
+        };
+        Self {
+            domain: ShortWeierstrass { a, b },
+            codomain,
+            s,
+        }
+    }
+
+    pub fn from_generator<C: Curve<F>>(curve: &C, generator: Point<F>, ell: u64) -> Self {
+        let iterations = ell / 2;
+        let mut kernel = Vec::with_capacity(iterations as usize);
+        let mut current = generator;
+        for _ in 0..iterations {
+            kernel.push(current);
+            current = curve.add(current, generator);
+        }
+        Self::velu(curve, kernel)
+    }
+
+    pub fn domain(&self) -> &ShortWeierstrass<F> {
+        &self.domain
+    }
+
+    pub fn codomain(&self) -> &ShortWeierstrass<F> {
+        &self.codomain
+    }
+
+    pub fn evaluate(&self, p: Point<F>) -> Point<F> {
+        let (x, y) = match p.representation() {
+            PointRepresentation::Infinity => return Point::infinity(),
+            PointRepresentation::Affine { x, y } => (x, y),
+        };
+        if self.s.iter().any(|(x_q, _, _)| *x_q == x) {
+            return Point::infinity();
+        }
+        let mut x_new = x;
+        let mut deriv_sum = FieldElement::zero();
+        for &(x_q, v_q, u_q) in &self.s {
+            let inv = (x - x_q).inv();
+            let inv2 = inv * inv;
+            let inv3 = inv2 * inv;
+            x_new += v_q * inv + u_q * inv2;
+            deriv_sum += v_q * inv2 + FieldElement::two() * u_q * inv3;
+        }
+        let y_new = y - y * deriv_sum;
+        Point::from_affine(x_new, y_new)
+    }
+}
+
+pub mod csidh {
+    use super::{Isogeny, ShortWeierstrass};
+    use curve::{Curve, Scalar};
+    use field::{Csidh512FieldOrder, Field, FieldElement};
+
+    pub const PRIMES: [u64; 74] = [
+        3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89,
+        97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181,
+        191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281,
+        283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 587,
+    ];
+
+    type Fp = Csidh512FieldOrder;
+    type Curve512 = ShortWeierstrass<Fp>;
+
+    pub fn base_curve() -> Curve512 {
+        ShortWeierstrass {
+            a: FieldElement::one(),
+            b: FieldElement::zero(),
+        }
+    }
+
+    fn twist(e: &Curve512) -> Curve512 {
+        ShortWeierstrass {
+            a: e.a,
+            b: FieldElement::zero() - e.b,
+        }
+    }
+
+    fn cofactor_for(ell: u64) -> Scalar<Fp> {
+        let mut limbs = [0u64; 8];
+        limbs[0] = 4;
+        for &l in PRIMES.iter() {
+            if l != ell {
+                mul_u64_in_place(&mut limbs, l);
+            }
+        }
+        Scalar::from_limbs(limbs)
+    }
+
+    fn mul_u64_in_place(limbs: &mut [u64; 8], factor: u64) {
+        let mut carry: u64 = 0;
+        for l in limbs.iter_mut() {
+            let (product, next_carry) = l.carrying_mul_add(factor, 0, carry);
+            *l = product;
+            carry = next_carry;
+        }
+        debug_assert_eq!(carry, 0, "cofactor overflowed 512 bits");
+    }
+
+    fn random_field<R: FnMut() -> u64>(next_u64: &mut R) -> FieldElement<Fp> {
+        loop {
+            let mut limbs = [0u64; 8];
+            for l in limbs.iter_mut() {
+                *l = next_u64();
+            }
+            if less_than_modulus(&limbs) {
+                return FieldElement::from_limbs_unchecked(limbs);
+            }
+        }
+    }
+
+    fn less_than_modulus(a: &[u64; 8]) -> bool {
+        let m = Fp::MODULUS;
+        for i in (0..8).rev() {
+            if a[i] != m[i] {
+                return a[i] < m[i];
+            }
+        }
+        false
+    }
+
+    pub fn action<R: FnMut() -> u64>(
+        start: &Curve512,
+        private_key: &[i8; PRIMES.len()],
+        next_u64: &mut R,
+    ) -> Curve512 {
+        let mut curve = *start;
+        let mut e = *private_key;
+
+        while e.iter().any(|&x| x != 0) {
+            for i in 0..PRIMES.len() {
+                if e[i] == 0 {
+                    continue;
+                }
+                let ell = PRIMES[i];
+                let sign: i8 = if e[i] > 0 { 1 } else { -1 };
+                let work = if sign > 0 { curve } else { twist(&curve) };
+                let cofactor = cofactor_for(ell);
+
+                let torsion = loop {
+                    let x = random_field(next_u64);
+                    let Some(p) = work.lift(x) else {
+                        continue;
+                    };
+                    let q = work.multiply(cofactor, p);
+                    if !q.is_infinity() {
+                        break q;
+                    }
+                };
+
+                let iso = Isogeny::from_generator(&work, torsion, ell);
+                let new_work = *iso.codomain();
+                curve = if sign > 0 {
+                    new_work
+                } else {
+                    twist(&new_work)
+                };
+                e[i] -= sign;
+            }
+        }
+        curve
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use curve::Secp256k1Curve;
+    use field::Secp256k1FieldOrder as Fp;
+
+    type Fe = FieldElement<Fp>;
+
+    #[test]
+    fn empty_kernel_yields_identity_isogeny() {
+        let iso = Isogeny::velu(&Secp256k1Curve, vec![]);
+        assert_eq!(iso.codomain().a, Secp256k1Curve::a());
+        assert_eq!(iso.codomain().b, Secp256k1Curve::b());
+        assert_eq!(
+            iso.evaluate(Secp256k1Curve::generator()),
+            Secp256k1Curve::generator()
+        );
+        assert!(iso.evaluate(Point::infinity()).is_infinity());
+    }
+
+    #[test]
+    fn empty_kernel_preserves_j_invariant() {
+        let iso = Isogeny::velu(&Secp256k1Curve, vec![]);
+        assert_eq!(iso.codomain().j_invariant(), Secp256k1Curve.j_invariant());
+    }
+
+    #[derive(Debug, Clone, Copy, PartialEq, Eq)]
+    struct E;
+    impl Curve<Fp> for E {
+        fn a(&self) -> Fe {
+            Fe::one()
+        }
+        fn b(&self) -> Fe {
+            Fe::zero() - Fe::two()
+        }
+    }
+
+    #[test]
+    fn two_isogeny_codomain_matches_hand_computation() {
+        let kernel = vec![Point::from_affine(Fe::one(), Fe::zero())];
+        let iso = Isogeny::velu(&E, kernel);
+        assert_eq!(iso.codomain().a, Fe::one() - Fe::from_u64(20));
+        assert_eq!(iso.codomain().b, (Fe::zero() - Fe::two()) - Fe::from_u64(28));
+    }
+
+    #[test]
+    fn two_isogeny_maps_kernel_to_infinity() {
+        let two_torsion = Point::from_affine(Fe::one(), Fe::zero());
+        let iso = Isogeny::velu(&E, vec![two_torsion]);
+        assert!(iso.evaluate(two_torsion).is_infinity());
+    }
+
+    #[test]
+    fn dedup_collapses_pm_pair_by_x_coordinate() {
+        let p = Point::from_affine(Fe::from_u64(42), Fe::from_u64(7));
+        let neg_p = Point::from_affine(Fe::from_u64(42), Fe::zero() - Fe::from_u64(7));
+        let iso = Isogeny::velu(&Secp256k1Curve, vec![p, neg_p]);
+        assert_eq!(iso.s.len(), 1);
+    }
+
+    #[test]
+    fn ignores_identity_in_kernel_list() {
+        let iso = Isogeny::velu(&Secp256k1Curve, vec![Point::infinity(), Point::infinity()]);
+        assert!(iso.s.is_empty());
+        assert_eq!(iso.codomain().a, Secp256k1Curve::a());
+        assert_eq!(iso.codomain().b, Secp256k1Curve::b());
+    }
+
+    #[test]
+    fn from_generator_two_isogeny_matches_velu() {
+        let two_torsion = Point::from_affine(Fe::one(), Fe::zero());
+        let iso_gen = Isogeny::from_generator(&E, two_torsion, 2);
+        let iso_velu = Isogeny::velu(&E, vec![two_torsion]);
+        assert_eq!(iso_gen.codomain().a, iso_velu.codomain().a);
+        assert_eq!(iso_gen.codomain().b, iso_velu.codomain().b);
+    }
+
+    #[test]
+    fn from_generator_odd_ell_enumerates_half() {
+        let g = Secp256k1Curve::generator();
+        assert_eq!(Isogeny::from_generator(&Secp256k1Curve, g, 3).s.len(), 1);
+        assert_eq!(Isogeny::from_generator(&Secp256k1Curve, g, 5).s.len(), 2);
+        assert_eq!(Isogeny::from_generator(&Secp256k1Curve, g, 7).s.len(), 3);
+    }
+}
+
+#[cfg(test)]
+mod csidh_tests {
+    use super::csidh::*;
+    use curve::Curve;
+    use field::{Csidh512FieldOrder, FieldElement};
+
+    #[test]
+    fn base_curve_is_y2_eq_x3_plus_x() {
+        let e0 = base_curve();
+        assert_eq!(e0.a, FieldElement::<Csidh512FieldOrder>::one());
+        assert_eq!(e0.b, FieldElement::<Csidh512FieldOrder>::zero());
+    }
+
+    #[test]
+    fn base_curve_j_invariant_is_1728() {
+        assert_eq!(
+            base_curve().j_invariant(),
+            FieldElement::<Csidh512FieldOrder>::from_u64(1728)
+        );
+    }
+
+    #[test]
+    fn action_with_zero_key_is_identity() {
+        let e0 = base_curve();
+        let key = [0i8; PRIMES.len()];
+        let result = action(&e0, &key, &mut || panic!("no RNG needed for zero key"));
+        assert_eq!(result.a, e0.a);
+        assert_eq!(result.b, e0.b);
+    }
+
+    #[test]
+    fn primes_count_is_74() {
+        assert_eq!(PRIMES.len(), 74);
+        assert_eq!(PRIMES[0], 3);
+        assert_eq!(PRIMES[PRIMES.len() - 1], 587);
+    }
+
+    #[test]
+    #[ignore = "slow — run with cargo test --release -- --ignored"]
+    fn action_commutes() {
+        use rand::{RngExt, SeedableRng, rngs::StdRng};
+
+        let mut rng = StdRng::seed_from_u64(0xC51D_DEAD_BEEF_0001);
+        let e0 = base_curve();
+        let mut a_key = [0i8; PRIMES.len()];
+        a_key[0] = 1;
+        let mut b_key = [0i8; PRIMES.len()];
+        b_key[1] = 1;
+
+        let a_pub = action(&e0, &a_key, &mut || rng.random::<u64>());
+        let b_pub = action(&e0, &b_key, &mut || rng.random::<u64>());
+        let a_shared = action(&b_pub, &a_key, &mut || rng.random::<u64>());
+        let b_shared = action(&a_pub, &b_key, &mut || rng.random::<u64>());
+
+        assert_eq!(a_shared.j_invariant(), b_shared.j_invariant());
+    }
+}