Repository: SpectralSequences/sseq (language: python)

Issue:
Support p = 11
We want to resolve modules related to Smith-Toda complexes, which are most relevant at primes bigger than 7.


Discussion / hints:

Fix patch:
diff --git a/ext/crates/algebra/benches/milnor.rs b/ext/crates/algebra/benches/milnor.rs
index ac7cf52fb..7025870c6 100644
--- a/ext/crates/algebra/benches/milnor.rs
+++ b/ext/crates/algebra/benches/milnor.rs
@@ -1,6 +1,6 @@
 use algebra::milnor_algebra::{PPartAllocation, PPartEntry, PPartMultiplier};
 use bencher::{benchmark_group, benchmark_main, Bencher};
-use fp::prime::ValidPrime;
+use fp::prime::{Prime, ValidPrime};
 
 fn ppart_inner<const MOD4: bool>(
     bench: &mut Bencher,
@@ -24,7 +24,7 @@ fn ppart_inner<const MOD4: bool>(
             if MOD4 {
                 assert!(c < 4);
             } else {
-                assert!(c < *p);
+                assert!(c < p.as_u32());
             }
         }
     });
diff --git a/ext/crates/algebra/src/algebra/adem_algebra.rs b/ext/crates/algebra/src/algebra/adem_algebra.rs
index 623ae4f2b..223f5c72c 100644
--- a/ext/crates/algebra/src/algebra/adem_algebra.rs
+++ b/ext/crates/algebra/src/algebra/adem_algebra.rs
@@ -5,7 +5,10 @@ use std::fmt;
 use itertools::Itertools;
 use rustc_hash::FxHashMap as HashMap;
 
-use fp::prime::{BinomialIterator, BitflagIterator, ValidPrime};
+use fp::prime::{
+    iter::{BinomialIterator, BitflagIterator},
+    Prime, ValidPrime,
+};
 use fp::vector::{FpVector, SliceMut};
 use once::OnceVec;
 
@@ -83,11 +86,10 @@ impl AdemBasisElement {
     fn excess(&self, p: ValidPrime) -> u32 {
         if self.degree < 2 {
             self.degree as u32 // special case these so that self.ps[0] is defined
-        } else if *p == 2 {
+        } else if p == 2 {
             2 * self.ps[0] - self.ps.iter().sum::<u32>()
         } else {
-            2 * (*p * self.ps[0] - (*p - 1) * self.ps.iter().sum::<u32>())
-                + 2 * (self.bocksteins & 1)
+            2 * (p * self.ps[0] - (p - 1) * self.ps.iter().sum::<u32>()) + 2 * (self.bocksteins & 1)
                 - self.bocksteins.count_ones()
         }
     }
@@ -181,9 +183,9 @@ impl Algebra for AdemAlgebra {
 
     fn magic(&self) -> u32 {
         if self.unstable_enabled {
-            1 + (*self.prime() << 16)
+            1 + (self.prime() << 16)
         } else {
-            *self.prime() << 16
+            self.prime() << 16
         }
     }
 
@@ -200,19 +202,19 @@ impl Algebra for AdemAlgebra {
                     degree: 1,
                     bocksteins: 1,
                     ps: vec![],
-                    p_or_sq: *self.prime() != 2,
+                    p_or_sq: self.prime() != 2,
                 },
             ));
             products.push((
                 "h_0".to_string(),
                 AdemBasisElement {
-                    degree: (2 * (*self.prime()) - 2) as i32,
+                    degree: (2 * self.prime() - 2) as i32,
                     bocksteins: 0,
                     ps: vec![1],
-                    p_or_sq: *self.prime() != 2,
+                    p_or_sq: self.prime() != 2,
                 },
             ));
-            (2 * (*self.prime()) - 2) as i32
+            (2 * self.prime() - 2) as i32
         } else {
             for i in 0..4 {
                 let degree = 1 << i; // degree is 2^hi
@@ -223,7 +225,7 @@ impl Algebra for AdemAlgebra {
                         degree,
                         bocksteins: 0,
                         ps,
-                        p_or_sq: *self.prime() != 2,
+                        p_or_sq: self.prime() != 2,
                     },
                 ));
             }
@@ -333,7 +335,7 @@ impl GeneratedAlgebra for AdemAlgebra {
             if degree == 1 {
                 "b".to_string()
             } else {
-                format!("P{}", degree as u32 / (2 * (*self.prime()) - 2))
+                format!("P{}", degree as u32 / (2 * self.prime() - 2))
             }
         } else {
             format!("Sq{degree}")
@@ -341,7 +343,7 @@ impl GeneratedAlgebra for AdemAlgebra {
     }
 
     fn generators(&self, degree: i32) -> Vec<usize> {
-        let p = *self.prime();
+        let p = self.prime();
         if degree == 0 {
             return vec![];
         }
@@ -365,7 +367,7 @@ impl GeneratedAlgebra for AdemAlgebra {
                 degree,
                 bocksteins: 0,
                 ps: vec![degree as u32 / (2 * p - 2)],
-                p_or_sq: *self.prime() != 2,
+                p_or_sq: self.prime() != 2,
             });
             vec![idx]
         } else {
@@ -377,7 +379,7 @@ impl GeneratedAlgebra for AdemAlgebra {
                 degree,
                 bocksteins: 0,
                 ps: vec![degree as u32],
-                p_or_sq: *self.prime() != 2,
+                p_or_sq: self.prime() != 2,
             });
             vec![idx]
         }
@@ -412,14 +414,14 @@ impl GeneratedAlgebra for AdemAlgebra {
             // Adem relation
             let first_sq = self.beps_pn(0, x);
             let second_sq = self.beps_pn(b, y);
-            relation.push((*p - 1, first_sq, second_sq));
+            relation.push((p - 1, first_sq, second_sq));
             for e1 in 0..=b {
                 let e2 = b - e1;
                 // e1 and e2 determine where a bockstein shows up.
                 // e1 determines if a bockstein shows up in front
                 // e2 determines if a bockstein shows up in middle
                 // So our output term looks like b^{e1} P^{x+y-j} b^{e2} P^{j}
-                for j in 0..=x / *p {
+                for j in 0..=x / p {
                     let c = combinatorics::adem_relation_coefficient(p, x, y, j, e1, e2);
                     if c == 0 {
                         continue;
@@ -432,7 +434,7 @@ impl GeneratedAlgebra for AdemAlgebra {
                             vec![x + y - j, j]
                         },
                         bocksteins: e1 | (e2 << 1),
-                        p_or_sq: *self.prime() != 2,
+                        p_or_sq: self.prime() != 2,
                     });
                     relation.push((c, (degree, idx), (0, 0)));
                 }
@@ -487,7 +489,7 @@ impl AdemAlgebra {
         let excess_table = OnceVec::new();
         Self {
             p,
-            generic: *p != 2,
+            generic: p != 2,
             unstable_enabled,
             even_basis_table,
             basis_table,
@@ -503,7 +505,7 @@ impl AdemAlgebra {
 
     pub fn q(&self) -> i32 {
         if self.generic {
-            2 * (*self.prime() as i32 - 1)
+            2 * (self.prime().as_i32() - 1)
         } else {
             1
         }
@@ -517,7 +519,7 @@ impl AdemAlgebra {
                 degree: 0,
                 bocksteins: if self.generic { u32::MAX << 1 } else { 0 },
                 ps: vec![],
-                p_or_sq: *p != 2,
+                p_or_sq: p != 2,
             }]
         });
 
@@ -530,7 +532,7 @@ impl AdemAlgebra {
                 degree: n,
                 bocksteins: if self.generic { u32::MAX << 2 } else { 0 },
                 ps: vec![n as u32],
-                p_or_sq: *self.prime() != 2,
+                p_or_sq: self.prime() != 2,
             });
 
             // last = last term. We append (last,) to the end of
@@ -540,12 +542,12 @@ impl AdemAlgebra {
             // or (p+1) * last <= n or last <= n/(p+1). We order the squares in decreasing
             // order of their last element so that as we walk over the previous basis
             // when we find a square whose end is too small, we can break.
-            for last in (1..=n / (*p as i32 + 1)).rev() {
+            for last in (1..=n / (p.as_i32() + 1)).rev() {
                 let previous_basis = &self.even_basis_table[(n - last) as usize];
                 for prev_elt in previous_basis {
                     let prev_elt_p_len = prev_elt.ps.len();
                     let old_last_sq = prev_elt.ps[prev_elt_p_len - 1] as i32;
-                    if old_last_sq < *p as i32 * last {
+                    if old_last_sq < p.as_i32() * last {
                         break;
                     }
                     // Write new basis element to basis element buffer
@@ -555,7 +557,7 @@ impl AdemAlgebra {
                     // A bit in bocksteins shall be set if it's illegal for a bockstein to occur there.
                     let mut bocksteins = prev_elt.bocksteins;
                     if self.generic {
-                        bocksteins |= if old_last_sq == *p as i32 * last {
+                        bocksteins |= if old_last_sq == p.as_i32() * last {
                             1 << prev_elt_p_len
                         } else {
                             0
@@ -569,7 +571,7 @@ impl AdemAlgebra {
                         degree,
                         bocksteins,
                         ps,
-                        p_or_sq: *p != 2,
+                        p_or_sq: p != 2,
                     });
                 }
             }
@@ -593,7 +595,7 @@ impl AdemAlgebra {
         self.generate_basis_even(max_degree);
 
         let p = self.prime();
-        let q = 2 * (*p as i32 - 1);
+        let q = 2 * (p.as_i32() - 1);
         self.basis_table.extend(max_degree as usize, |n| {
             let n = n as i32;
             let residue = n % q;
@@ -626,7 +628,7 @@ impl AdemAlgebra {
                             degree: n,
                             bocksteins,
                             ps: elt.ps.clone(),
-                            p_or_sq: *p != 2,
+                            p_or_sq: p != 2,
                         });
                         // BinomialIterator does not have a stopping condition, so we need to stop
                         // manually here
@@ -778,7 +780,7 @@ impl AdemAlgebra {
     fn generate_multiplication_table_generic(&self, max_degree: i32) {
         self.multiplication_table.extend(0, |_| Vec::new());
 
-        let q = 2 * (*self.prime()) as i32 - 2;
+        let q = 2 * self.prime().as_i32() - 2;
         self.multiplication_table.extend(max_degree as usize, |n| {
             let n = n as i32;
             let mut table: Vec<Vec<FpVector>> = Vec::with_capacity(2 * (n / q + 1) as usize);
@@ -821,7 +823,7 @@ impl AdemAlgebra {
         x: i32,
         idx: usize,
     ) -> FpVector {
-        let p: i32 = *self.prime() as i32; // we use p for the i32 version and self.p for the u32 version
+        let p: i32 = self.prime().as_i32(); // we use p for the i32 version and self.p for the u32 version
         let q: i32 = 2 * p - 2;
 
         let x: u32 = x as u32;
@@ -864,7 +866,7 @@ impl AdemAlgebra {
             .clone();
 
         let b: u32 = working_elt.bocksteins & 1;
-        if working_elt.ps.is_empty() || i >= (*self.prime()) * working_elt.ps[0] + b {
+        if working_elt.ps.is_empty() || i >= self.prime() * working_elt.ps[0] + b {
             working_elt.ps.insert(0, i);
             working_elt.bocksteins <<= 1;
             working_elt.degree = n;
@@ -881,7 +883,7 @@ impl AdemAlgebra {
 
         if b == 0 {
             // We use P^i P^j = \sum ... P^{i + j - k} P^k
-            for k in 0..=i / (*self.prime()) {
+            for k in 0..=i / self.prime() {
                 let c = combinatorics::adem_relation_coefficient(self.prime(), i, j, k, 0, 0);
                 if c == 0 {
                     continue;
@@ -899,7 +901,7 @@ impl AdemAlgebra {
                     [(n - q * (i + j - k) as i32) as usize][2 * k as usize][tail_idx];
                 for (id, coeff) in rest_reduced.iter().enumerate() {
                     let source = &table[2 * (i + j - k) as usize][id];
-                    result.add(source, (c * coeff) % *self.prime());
+                    result.add(source, (c * coeff) % self.prime());
                 }
             }
         } else {
@@ -919,7 +921,7 @@ impl AdemAlgebra {
                 result.add_basis_element(index, c);
             }
 
-            for k in 1..=i / (*self.prime()) {
+            for k in 1..=i / self.prime() {
                 // \beta P^{i + j - k} P^k
                 let c = combinatorics::adem_relation_coefficient(self.prime(), i, j, k, 1, 0);
                 if c != 0 {
@@ -927,7 +929,7 @@ impl AdemAlgebra {
                         [(n - q * (i + j - k) as i32 - 1) as usize][2 * k as usize][tail_idx];
                     for (id, coeff) in rest_reduced.iter().enumerate() {
                         let source = &table[1 + 2 * (i + j - k) as usize][id];
-                        result.add(source, (c * coeff) % *self.prime());
+                        result.add(source, (c * coeff) % self.prime());
                     }
                 }
 
@@ -938,7 +940,7 @@ impl AdemAlgebra {
                         [(n - q * (i + j - k) as i32) as usize][1 + 2 * k as usize][tail_idx];
                     for (id, coeff) in rest_reduced.iter().enumerate() {
                         let source = &table[2 * (i + j - k) as usize][id];
-                        result.add(source, (c * coeff) % *self.prime());
+                        result.add(source, (c * coeff) % self.prime());
                     }
                 }
             }
@@ -975,7 +977,7 @@ impl AdemAlgebra {
             degree: r.degree + s.degree,
             bocksteins: 0,
             ps: Vec::with_capacity(r.ps.len() + s.ps.len()),
-            p_or_sq: *self.prime() != 2,
+            p_or_sq: self.prime() != 2,
         };
         if self.generic && (r.bocksteins >> r.ps.len()) & s.bocksteins & 1 == 1 {
             // If there is a bockstein at the end of r and one at the beginning of s, these run into each other
@@ -1026,7 +1028,7 @@ impl AdemAlgebra {
         excess: i32,
     ) {
         let q = if self.generic {
-            2 * (*self.prime()) - 2
+            2 * self.prime() - 2
         } else {
             1
         };
@@ -1095,7 +1097,7 @@ impl AdemAlgebra {
             degree: monomial.degree,
             bocksteins: 0,
             ps: monomial.ps[0..idx].to_vec(),
-            p_or_sq: *self.prime() != 2,
+            p_or_sq: self.prime() != 2,
         };
 
         for (it_idx, _value) in reduced_tail.iter_nonzero() {
@@ -1123,7 +1125,7 @@ impl AdemAlgebra {
         excess: i32,
         stop_early: bool,
     ) {
-        let p = *self.prime();
+        let p = self.prime();
         let q = 2 * p - 2;
         // Check for admissibility
         let b1 = if idx >= 0 {
@@ -1160,7 +1162,7 @@ impl AdemAlgebra {
             degree: monomial.degree,
             bocksteins: 0,
             ps: monomial.ps[0..idx].to_vec(),
-            p_or_sq: *self.prime() != 2,
+            p_or_sq: self.prime() != 2,
         };
 
         for (it_idx, it_value) in reduced_tail.iter_nonzero() {
@@ -1196,13 +1198,13 @@ impl AdemAlgebra {
                 degree: degree_first,
                 bocksteins: 0,
                 ps: vec![b.ps[0]],
-                p_or_sq: *self.prime() != 2,
+                p_or_sq: self.prime() != 2,
             });
             let idx_rest = self.basis_element_to_index(&AdemBasisElement {
                 degree: degree_rest,
                 bocksteins: 0,
                 ps: ps_rest,
-                p_or_sq: *self.prime() != 2,
+                p_or_sq: self.prime() != 2,
             });
             return vec![(1, (degree_first, idx_first), (degree_rest, idx_rest))];
         }
@@ -1216,13 +1218,13 @@ impl AdemAlgebra {
             degree: first_degree,
             bocksteins: 0,
             ps: vec![first_sq],
-            p_or_sq: *self.prime() != 2,
+            p_or_sq: self.prime() != 2,
         });
         let second_idx = self.basis_element_to_index(&AdemBasisElement {
             degree: second_degree,
             bocksteins: 0,
             ps: vec![second_sq],
-            p_or_sq: *self.prime() != 2,
+            p_or_sq: self.prime() != 2,
         });
         let mut out_vec = FpVector::new(fp::prime::TWO, self.dimension(degree));
         self.multiply_basis_elements(
@@ -1260,7 +1262,7 @@ impl AdemAlgebra {
             return vec![(1, (first_degree, first_idx), (rest_degree, rest_idx))];
         }
         if b.bocksteins != 0 || b.ps.len() != 1 {
-            let first_degree = (b.ps[0] * 2 * (*p - 1)) as i32;
+            let first_degree = (b.ps[0] * 2 * (p - 1)) as i32;
             let rest_degree = b.degree - first_degree;
             let ps_first = vec![b.ps[0]];
             let ps_rest = b.ps[1..].to_vec();
@@ -1268,13 +1270,13 @@ impl AdemAlgebra {
                 degree: first_degree,
                 bocksteins: 0,
                 ps: ps_first,
-                p_or_sq: *self.prime() != 2,
+                p_or_sq: self.prime() != 2,
             };
             let rest = AdemBasisElement {
                 degree: rest_degree,
                 bocksteins: b.bocksteins >> 1,
                 ps: ps_rest,
-                p_or_sq: *self.prime() != 2,
+                p_or_sq: self.prime() != 2,
             };
             let first_idx = self.basis_element_to_index(&first);
             let rest_idx = self.basis_element_to_index(&rest);
@@ -1285,27 +1287,27 @@ impl AdemAlgebra {
         let mut pow = 1;
         {
             let mut temp_sq = sq;
-            while temp_sq % *p == 0 {
-                temp_sq /= *p;
-                pow *= *p;
+            while temp_sq % p == 0 {
+                temp_sq /= p;
+                pow *= p;
             }
         }
 
         let first_sq = pow;
         let second_sq = sq - first_sq;
-        let first_degree = (first_sq * 2 * (*p - 1)) as i32;
-        let second_degree = (second_sq * 2 * (*p - 1)) as i32;
+        let first_degree = (first_sq * 2 * (p - 1)) as i32;
+        let second_degree = (second_sq * 2 * (p - 1)) as i32;
         let first_idx = self.basis_element_to_index(&AdemBasisElement {
             degree: first_degree,
             bocksteins: 0,
             ps: vec![first_sq],
-            p_or_sq: *self.prime() != 2,
+            p_or_sq: self.prime() != 2,
         });
         let second_idx = self.basis_element_to_index(&AdemBasisElement {
             degree: second_degree,
             bocksteins: 0,
             ps: vec![second_sq],
-            p_or_sq: *self.prime() != 2,
+            p_or_sq: self.prime() != 2,
         });
         let mut out_vec = FpVector::new(p, self.dimension(degree));
         self.multiply_basis_elements(
@@ -1319,16 +1321,16 @@ impl AdemAlgebra {
         let mut result = Vec::new();
         let c = out_vec.entry(idx);
         assert!(c != 0);
-        let c_inv = fp::prime::inverse(p, *p - c);
+        let c_inv = fp::prime::inverse(p, p - c);
         result.push((
-            ((*p - 1) * c_inv) % *p,
+            ((p - 1) * c_inv) % p,
             (first_degree, first_idx),
             (second_degree, second_idx),
         ));
         out_vec.set_entry(idx, 0);
         for (i, v) in out_vec.iter_nonzero() {
             let (c, t1, t2) = self.decompose_basis_element_generic(degree, i)[0];
-            result.push(((c_inv * c * v) % *p, t1, t2));
+            result.push(((c_inv * c * v) % p, t1, t2));
         }
         result
     }
@@ -1340,14 +1342,14 @@ impl AdemAlgebra {
             return (0, 0);
         }
 
-        let p = *self.prime();
+        let p = self.prime();
         let q = if self.generic { 2 * p - 2 } else { 1 };
         let degree = (x * q + e) as i32;
         let index = self.basis_element_to_index(&AdemBasisElement {
             degree,
             bocksteins: e,
             ps: vec![x],
-            p_or_sq: *self.prime() != 2,
+            p_or_sq: self.prime() != 2,
         });
         (degree, index)
     }
@@ -1407,7 +1409,7 @@ impl Bialgebra for AdemAlgebra {
             if op_deg == 1 {
                 vec![(1, 0, 0, 0), (0, 0, 1, 0)]
             } else {
-                let q = *self.prime() * 2 - 2;
+                let q = self.prime() * 2 - 2;
                 let op_deg = op_deg as u32;
                 assert_eq!(op_deg % q, 0);
 
diff --git a/ext/crates/algebra/src/algebra/algebra_trait.rs b/ext/crates/algebra/src/algebra/algebra_trait.rs
index e58a8a9e0..532b39167 100644
--- a/ext/crates/algebra/src/algebra/algebra_trait.rs
+++ b/ext/crates/algebra/src/algebra/algebra_trait.rs
@@ -83,7 +83,7 @@ pub trait Algebra: std::fmt::Display + Send + Sync + 'static {
         for (i, v) in s.iter_nonzero() {
             self.multiply_basis_elements(
                 result.copy(),
-                (coeff * v) % *p,
+                (coeff * v) % p,
                 r_degree,
                 r_idx,
                 s_degree,
@@ -109,7 +109,7 @@ pub trait Algebra: std::fmt::Display + Send + Sync + 'static {
         for (i, v) in r.iter_nonzero() {
             self.multiply_basis_elements(
                 result.copy(),
-                (coeff * v) % *p,
+                (coeff * v) % p,
                 r_degree,
                 i,
                 s_degree,
@@ -135,7 +135,7 @@ pub trait Algebra: std::fmt::Display + Send + Sync + 'static {
         for (i, v) in s.iter_nonzero() {
             self.multiply_element_by_basis_element(
                 result.copy(),
-                (coeff * v) % *p,
+                (coeff * v) % p,
                 r_degree,
                 r,
                 s_degree,
@@ -222,7 +222,7 @@ pub trait UnstableAlgebra: Algebra {
         for (i, v) in s.iter_nonzero() {
             self.multiply_basis_elements_unstable(
                 result.copy(),
-                (coeff * v) % *p,
+                (coeff * v) % p,
                 r_degree,
                 r_idx,
                 s_degree,
@@ -250,7 +250,7 @@ pub trait UnstableAlgebra: Algebra {
         for (i, v) in r.iter_nonzero() {
             self.multiply_basis_elements_unstable(
                 result.copy(),
-                (coeff * v) % *p,
+                (coeff * v) % p,
                 r_degree,
                 i,
                 s_degree,
@@ -278,7 +278,7 @@ pub trait UnstableAlgebra: Algebra {
         for (i, v) in s.iter_nonzero() {
             self.multiply_element_by_basis_element_unstable(
                 result.copy(),
-                (coeff * v) % *p,
+                (coeff * v) % p,
                 r_degree,
                 r,
                 s_degree,
diff --git a/ext/crates/algebra/src/algebra/combinatorics.rs b/ext/crates/algebra/src/algebra/combinatorics.rs
index 79c8661dc..5983d0a70 100644
--- a/ext/crates/algebra/src/algebra/combinatorics.rs
+++ b/ext/crates/algebra/src/algebra/combinatorics.rs
@@ -1,7 +1,7 @@
 use fp::vector::FpVector;
 use once::OnceVec;
 
-use fp::prime::{minus_one_to_the_n, Binomial, ValidPrime};
+use fp::prime::{minus_one_to_the_n, Binomial, Prime, ValidPrime};
 use fp::{const_for, MAX_MULTINOMIAL_LEN, NUM_PRIMES, PRIMES, PRIME_TO_INDEX_MAP};
 
 pub const MAX_XI_TAU: usize = MAX_MULTINOMIAL_LEN;
@@ -39,7 +39,7 @@ const TAU_DEGREES: [[i32; MAX_XI_TAU]; NUM_PRIMES] = {
 };
 
 pub fn adem_relation_coefficient(p: ValidPrime, x: u32, y: u32, j: u32, e1: u32, e2: u32) -> u32 {
-    let pi32 = *p as i32;
+    let pi32 = p.as_i32();
     let x = x as i32;
     let y = y as i32;
     let j = j as i32;
@@ -49,12 +49,11 @@ pub fn adem_relation_coefficient(p: ValidPrime, x: u32, y: u32, j: u32, e1: u32,
     if c == 0 {
         return 0;
     }
-    c *= minus_one_to_the_n(*p, (x + j) + e2);
-    c % *p
+    c *= minus_one_to_the_n(p, (x + j) + e2);
+    c % p
 }
 
 pub fn inadmissible_pairs(p: ValidPrime, generic: bool, degree: i32) -> Vec<(u32, u32, u32)> {
-    let p = *p;
     let degree = degree as u32;
     let q = if generic { 2 * p - 2 } else { 1 };
     // (i, b, j) means P^i P^j if b = 0, or P^i b P^j if b = 1.
@@ -83,11 +82,11 @@ pub fn inadmissible_pairs(p: ValidPrime, generic: bool, degree: i32) -> Vec<(u32
 }
 
 pub fn tau_degrees(p: ValidPrime) -> &'static [i32] {
-    &TAU_DEGREES[PRIME_TO_INDEX_MAP[*p as usize]]
+    &TAU_DEGREES[PRIME_TO_INDEX_MAP[p.as_usize()]]
 }
 
 pub fn xi_degrees(p: ValidPrime) -> &'static [i32] {
-    &XI_DEGREES[PRIME_TO_INDEX_MAP[*p as usize]]
+    &XI_DEGREES[PRIME_TO_INDEX_MAP[p.as_usize()]]
 }
 
 pub struct TruncatedPolynomialMonomialBasis {
@@ -142,7 +141,7 @@ impl TruncatedPolynomialMonomialBasis {
 
     pub fn add_gens_and_calculate_parts(&self, degree: i32, new_gens: usize) {
         assert!(degree as usize == self.gens.len());
-        let p = *self.p;
+        let p = self.p;
         let idx = self.gens[degree as usize - 1].0 + self.gens[degree as usize - 1].1;
         self.gens.push((idx, new_gens));
         let mut new_parts_by_max = Vec::new();
diff --git a/ext/crates/algebra/src/algebra/milnor_algebra.rs b/ext/crates/algebra/src/algebra/milnor_algebra.rs
index 08aa2a101..aaad79eb9 100644
--- a/ext/crates/algebra/src/algebra/milnor_algebra.rs
+++ b/ext/crates/algebra/src/algebra/milnor_algebra.rs
@@ -1,10 +1,11 @@
+use fp::prime::Prime;
 use itertools::Itertools;
 use rustc_hash::FxHashMap as HashMap;
 use std::cell::Cell;
 
 use crate::algebra::combinatorics;
 use crate::algebra::{Algebra, Bialgebra, GeneratedAlgebra, UnstableAlgebra};
-use fp::prime::{factor_pk, integer_power, Binomial, BitflagIterator, ValidPrime};
+use fp::prime::{factor_pk, iter::BitflagIterator, Binomial, ValidPrime};
 use fp::vector::{FpVector, Slice, SliceMut};
 use once::OnceVec;
 
@@ -130,7 +131,7 @@ impl MilnorBasisElement {
     }
 
     fn excess(&self, p: ValidPrime) -> u32 {
-        if *p == 2 {
+        if p == 2 {
             self.p_part.iter().sum::<PPartEntry>()
         } else {
             self.q_part.count_ones() + 2 * self.p_part.iter().sum::<PPartEntry>()
@@ -146,7 +147,7 @@ impl MilnorBasisElement {
 
     /// Update the degree component to the correct degree
     pub fn compute_degree(&mut self, p: ValidPrime) {
-        let q = if *p == 2 { 1 } else { 2 * (*p as i32 - 1) };
+        let q = if p == 2 { 1 } else { 2 * (p.as_i32() - 1) };
         let xi_degrees = combinatorics::xi_degrees(p);
         let tau_degrees = combinatorics::tau_degrees(p);
 
@@ -298,7 +299,7 @@ impl MilnorAlgebra {
         Self {
             p,
             #[cfg(feature = "odd-primes")]
-            generic: *p != 2,
+            generic: p != 2,
             unstable_enabled,
             profile,
             ppart_table: OnceVec::new(),
@@ -325,7 +326,7 @@ impl MilnorAlgebra {
 
     pub fn q(&self) -> i32 {
         if self.generic() {
-            2 * (*self.prime() as i32 - 1)
+            2 * (self.prime().as_i32() - 1)
         } else {
             1
         }
@@ -362,7 +363,7 @@ impl Algebra for MilnorAlgebra {
     }
 
     fn magic(&self) -> u32 {
-        (*self.p << 16)
+        (self.p << 16)
             + if self.profile.is_trivial() {
                 0x8000
             } else {
@@ -394,13 +395,13 @@ impl Algebra for MilnorAlgebra {
                 products.push((
                     "h_0".to_string(),
                     MilnorBasisElement {
-                        degree: (2 * (*self.prime()) - 2) as i32,
+                        degree: (2 * self.prime() - 2) as i32,
                         q_part: 0,
                         p_part: vec![1],
                     },
                 ));
             }
-            max_degree = (2 * (*self.prime()) - 2) as i32;
+            max_degree = (2 * self.prime() - 2) as i32;
         } else {
             let mut max = 4;
             if !self.profile.p_part.is_empty() {
@@ -540,7 +541,7 @@ impl Algebra for MilnorAlgebra {
             for (i, v) in s.iter_nonzero() {
                 allocation = self.multiply_with_allocation(
                     result.copy(),
-                    (coeff * v) % *p,
+                    (coeff * v) % p,
                     r,
                     self.basis_element_from_index(s_degree, i),
                     i32::MAX,
@@ -712,7 +713,7 @@ impl GeneratedAlgebra for MilnorAlgebra {
             return vec![0]; // Q_0
         }
 
-        let p = *self.prime();
+        let p = self.prime();
 
         // Check for the Q_k
         if self.generic() && degree % 2 == 1 {
@@ -806,7 +807,7 @@ impl GeneratedAlgebra for MilnorAlgebra {
                 let (first_degree, first_index) = self.try_beps_pn(0, x as PPartEntry)?;
                 let (second_degree, second_index) = self.try_beps_pn(b, y as PPartEntry)?;
                 relation.push((
-                    *p - 1,
+                    p - 1,
                     (first_degree, first_index),
                     (second_degree, second_index),
                 ));
@@ -816,7 +817,7 @@ impl GeneratedAlgebra for MilnorAlgebra {
                     // e1 determines whether a bockstein shows up in front
                     // e2 determines whether a bockstein shows up in middle
                     // So our output term looks like b^{e1} P^{x+y-j} b^{e2} P^{j}
-                    for j in 0..=x / *p {
+                    for j in 0..=x / p {
                         let c = combinatorics::adem_relation_coefficient(p, x, y, j, e1, e2);
                         if c == 0 {
                             continue;
@@ -847,7 +848,7 @@ impl MilnorAlgebra {
     fn compute_ppart(&self, max_degree: i32) {
         self.ppart_table.extend(0, |_| vec![Vec::new()]);
 
-        let p = *self.prime() as i32;
+        let p = self.prime().as_i32();
         let q = if p == 2 { 1 } else { 2 * p - 2 };
         let new_deg = max_degree / q;
 
@@ -855,8 +856,7 @@ impl MilnorAlgebra {
         let mut profile_list = Vec::with_capacity(xi_degrees.len());
         for i in 0..xi_degrees.len() {
             if i < self.profile.p_part.len() {
-                profile_list
-                    .push((integer_power(*self.prime(), self.profile.p_part[i]) - 1) as PPartEntry);
+                profile_list.push((self.prime().pow(self.profile.p_part[i]) - 1) as PPartEntry);
             } else if self.profile.truncated {
                 profile_list.push(0);
             } else {
@@ -898,7 +898,7 @@ impl MilnorAlgebra {
     }
 
     fn generate_basis_generic(&self, max_degree: i32) {
-        let q = 2 * (*self.prime()) - 2;
+        let q = 2 * self.prime() - 2;
         let tau_degrees = combinatorics::tau_degrees(self.prime());
 
         self.basis_table.extend(max_degree as usize, |d| {
@@ -1004,7 +1004,7 @@ impl MilnorAlgebra {
 
         for k in BitflagIterator::set_bit_iterator(f as u64) {
             let k = k as u32;
-            let pk = integer_power(*self.p, k) as PPartEntry;
+            let pk = self.p.pow(k) as PPartEntry;
             std::mem::swap(&mut new_result, &mut old_result);
             new_result.clear();
 
@@ -1048,7 +1048,7 @@ impl MilnorAlgebra {
                     let c = if larger_q % 2 == 0 {
                         *coef
                     } else {
-                        *coef * (*self.prime() - 1)
+                        *coef * (self.prime() - 1)
                     };
 
                     new_result.push((c, m));
@@ -1271,7 +1271,7 @@ impl<'a, const MOD4: bool> PPartMultiplier<'a, MOD4> {
         degree: i32,
     ) -> Self {
         if MOD4 {
-            assert_eq!(*p, 2);
+            assert_eq!(p, 2);
         }
         let rows = r.len() + 1;
         let cols = s.len() + 1;
@@ -1327,7 +1327,7 @@ impl<'a, const MOD4: bool> PPartMultiplier<'a, MOD4> {
     ///
     /// TODO: Improve odd prime performance
     fn next_val(&self, sum: PPartEntry, k: PPartEntry, max: PPartEntry) -> PPartEntry {
-        match *self.prime() {
+        match self.prime().as_u32() {
             2 => {
                 if MOD4 {
                     // x.count_ones() + y.count_ones() - (x + y).count_ones() is the number of
@@ -1369,7 +1369,7 @@ impl<'a, const MOD4: bool> PPartMultiplier<'a, MOD4> {
             let mut total = self.M[i][0];
             let mut p_to_the_j = 1;
             for j in 1..self.cols {
-                p_to_the_j *= *self.prime() as PPartEntry;
+                p_to_the_j *= self.prime().as_u32() as PPartEntry;
                 if total < p_to_the_j {
                     // We don't have enough weight left in the entries above this one in the column to increment this cell.
                     // Add the weight from this cell to the total, we can use it to increment a cell lower down.
@@ -1430,7 +1430,7 @@ impl<'a, const MOD4: bool> Iterator for PPartMultiplier<'a, MOD4> {
     type Item = u32;
 
     fn next(&mut self) -> Option<u32> {
-        let p = *self.prime() as PPartEntry;
+        let p = self.prime().as_u32() as PPartEntry;
         'outer: loop {
             self.ans.p_part.clear();
             let mut coef = 1;
@@ -1483,7 +1483,7 @@ impl<'a, const MOD4: bool> Iterator for PPartMultiplier<'a, MOD4> {
                     let i_max = std::cmp::min(1 + diag_idx, self.rows);
                     let mut sum = 0;
 
-                    if *self.prime() == 2 {
+                    if self.prime() == 2 {
                         if MOD4 {
                             for i in i_min..i_max {
                                 let entry = self.M[i][diag_idx - i];
@@ -1550,10 +1550,10 @@ impl MilnorAlgebra {
         let i = basis.q_part.trailing_zeros();
         // If the basis element is just Q_{k+1}, we decompose Q_{k+1} = P(p^k) Q_k - Q_k P(p^k).
         if basis.q_part == 1 << i && basis.p_part.is_empty() {
-            let ppow = fp::prime::integer_power(*self.prime(), i - 1);
+            let ppow = self.prime().pow(i - 1);
 
             let q_degree = (2 * ppow - 1) as i32;
-            let p_degree = (ppow * (2 * (*self.prime()) - 2)) as i32;
+            let p_degree = (ppow * (2 * self.prime() - 2)) as i32;
 
             let p_idx = self
                 .basis_element_to_index(&MilnorBasisElement::from_p(
@@ -1572,7 +1572,7 @@ impl MilnorAlgebra {
 
             vec![
                 (1, (p_degree, p_idx), (q_degree, q_idx)),
-                (*self.prime() - 1, (q_degree, q_idx), (p_degree, p_idx)),
+                (self.prime() - 1, (q_degree, q_idx), (p_degree, p_idx)),
             ]
         } else {
             // Otherwise, separate out the first Q_k.
@@ -1616,7 +1616,7 @@ impl MilnorAlgebra {
 
         // out_vec is the remaining items we have to kill. We are done when this hits zero.
         let mut out_vec = FpVector::new(p, self.dimension(degree));
-        out_vec.set_entry(idx, *p - 1);
+        out_vec.set_entry(idx, p - 1);
 
         while let Some((idx, c)) = out_vec.iter_nonzero().next() {
             let b = self.basis_element_from_index(degree, idx);
@@ -1625,12 +1625,12 @@ impl MilnorAlgebra {
             if b.p_part[0..len - 1].iter().all(|&x| x == 0) {
                 // There is only one entry
                 let entry = b.p_part[len - 1];
-                let (k, m) = factor_pk(*p, entry);
+                let (k, m) = factor_pk(p, entry);
 
                 // This is a power of p
                 if m == 1 {
                     if len == 1 || !self.profile.is_an(self.generic()) {
-                        buffer.extend([(*p - c, (degree, idx), (0, 0))]);
+                        buffer.extend([(p - c, (degree, idx), (0, 0))]);
                     } else {
                         // Write this as [P(p^(len + k - 1)), P(0, .., 0, P^k)] plus higher order
                         // terms.
@@ -1655,7 +1655,7 @@ impl MilnorAlgebra {
                             degree: r_degree,
                         });
                         buffer.extend(vec![
-                            (*p - c, (l_degree, l_index), (r_degree, r_index)),
+                            (p - c, (l_degree, l_index), (r_degree, r_index)),
                             (c, (r_degree, r_index), (l_degree, l_index)),
                         ])
                     }
@@ -1681,7 +1681,7 @@ impl MilnorAlgebra {
                     let second = (elt.degree, self.basis_element_to_index(&elt));
 
                     let coef =
-                        *p - fp::prime::inverse(p, PPartEntry::binomial(p, pk + rem_entry, pk));
+                        p - fp::prime::inverse(p, PPartEntry::binomial(p, pk + rem_entry, pk));
                     buffer.extend([(coef, first, second)])
                 }
             } else {
@@ -1703,7 +1703,7 @@ impl MilnorAlgebra {
                     elt.p_part.pop();
                 }
                 let second = (elt.degree, self.basis_element_to_index(&elt));
-                buffer.extend([(*p - c, first, second)]);
+                buffer.extend([(p - c, first, second)]);
             };
             for (c, first, second) in &buffer {
                 self.multiply_basis_elements(
@@ -1738,7 +1738,7 @@ impl MilnorAlgebra {
 
 impl Bialgebra for MilnorAlgebra {
     fn coproduct(&self, op_deg: i32, op_idx: usize) -> Vec<(i32, usize, i32, usize)> {
-        assert_eq!(*self.prime(), 2, "Coproduct at odd primes not supported");
+        assert_eq!(self.prime(), 2, "Coproduct at odd primes not supported");
         if op_deg == 0 {
             return vec![(0, 0, 0, 0)];
         }
diff --git a/ext/crates/algebra/src/algebra/polynomial_algebra.rs b/ext/crates/algebra/src/algebra/polynomial_algebra.rs
index 89f855418..3dd362830 100644
--- a/ext/crates/algebra/src/algebra/polynomial_algebra.rs
+++ b/ext/crates/algebra/src/algebra/polynomial_algebra.rs
@@ -2,7 +2,7 @@ use itertools::Itertools;
 use rustc_hash::FxHashMap as HashMap;
 use std::fmt;
 
-use fp::prime::ValidPrime;
+use fp::prime::{Prime, ValidPrime};
 use fp::vector::{FpVector, SliceMut};
 use once::OnceVec;
 
@@ -123,7 +123,7 @@ pub trait PolynomialAlgebra: std::fmt::Display + Sized + Send + Sync + 'static {
     }
 
     fn frobenius_monomial(&self, target: &mut FpVector, source: &FpVector) {
-        let p = *self.prime() as i32;
+        let p = self.prime().as_i32();
         for (i, c) in source.iter_nonzero() {
             let (gen_degree, gen_index) = self.polynomial_monomials().internal_idx_to_gen_deg(i);
             let frob = self.frobenius_on_generator(gen_degree, gen_index);
@@ -141,7 +141,7 @@ pub trait PolynomialAlgebra: std::fmt::Display + Sized + Send + Sync + 'static {
         target: &mut PolynomialAlgebraMonomial,
         source: &PolynomialAlgebraMonomial,
     ) -> Option<u32> {
-        let minus_one = *self.prime() - 1;
+        let minus_one = self.prime() - 1;
         self.set_monomial_degree(target, target.degree + source.degree);
         let mut temp_source_ext = source.ext.clone();
         temp_source_ext.set_scratch_vector_size(target.ext.len());
@@ -177,7 +177,7 @@ pub trait PolynomialAlgebra: std::fmt::Display + Sized + Send + Sync + 'static {
         right_degree: i32,
         right: &FpVector,
     ) {
-        let p = *self.prime();
+        let p = self.prime();
         target.set_scratch_vector_size(self.dimension(left_degree + right_degree));
         for (left_idx, left_entry) in left.iter_nonzero() {
             for (right_idx, right_entry) in right.iter_nonzero() {
@@ -200,7 +200,7 @@ pub trait PolynomialAlgebra: std::fmt::Display + Sized + Send + Sync + 'static {
         left: &FpVector,
         right_mono: &PolynomialAlgebraMonomial,
     ) {
-        let p = *self.prime();
+        let p = self.prime();
         target.extend_len(self.dimension(left_degree + right_mono.degree));
         for (left_idx, left_entry) in left.iter_nonzero() {
             let mut target_mono = self.index_to_monomial(left_degree, left_idx).clone(); // Could reduce cloning a bit but probably best not to worry.
@@ -221,7 +221,7 @@ pub trait PolynomialAlgebra: std::fmt::Display + Sized + Send + Sync + 'static {
         right_degree: i32,
         right: &FpVector,
     ) {
-        let p = *self.prime();
+        let p = self.prime();
         target.extend_len(self.dimension(right_degree + left_mono.degree));
         for (right_idx, right_entry) in right.iter_nonzero() {
             let mut target_mono = left_mono.clone(); // Could reduce cloning a bit but probably best not to worry.
diff --git a/ext/crates/algebra/src/module/homomorphism/mod.rs b/ext/crates/algebra/src/module/homomorphism/mod.rs
index 733c5f532..ac7043832 100644
--- a/ext/crates/algebra/src/module/homomorphism/mod.rs
+++ b/ext/crates/algebra/src/module/homomorphism/mod.rs
@@ -71,7 +71,7 @@ pub trait ModuleHomomorphism: Send + Sync {
     fn apply(&self, mut result: SliceMut, coeff: u32, input_degree: i32, input: Slice) {
         let p = self.prime();
         for (i, v) in input.iter_nonzero() {
-            self.apply_to_basis_element(result.copy(), (coeff * v) % *p, input_degree, i);
+            self.apply_to_basis_element(result.copy(), (coeff * v) % p, input_degree, i);
         }
     }
 
diff --git a/ext/crates/algebra/src/module/module_trait.rs b/ext/crates/algebra/src/module/module_trait.rs
index 321dbe3e7..9947ac6d4 100644
--- a/ext/crates/algebra/src/module/module_trait.rs
+++ b/ext/crates/algebra/src/module/module_trait.rs
@@ -110,7 +110,7 @@ pub trait Module: std::fmt::Display + std::any::Any + Send + Sync {
         for (i, v) in input.iter_nonzero() {
             self.act_on_basis(
                 result.copy(),
-                (coeff * v) % *p,
+                (coeff * v) % p,
                 op_degree,
                 op_index,
                 input_degree,
@@ -134,7 +134,7 @@ pub trait Module: std::fmt::Display + std::any::Any + Send + Sync {
         for (i, v) in op.iter_nonzero() {
             self.act(
                 result.copy(),
-                (coeff * v) % *p,
+                (coeff * v) % p,
                 op_degree,
                 i,
                 input_degree,
@@ -157,7 +157,7 @@ pub trait Module: std::fmt::Display + std::any::Any + Send + Sync {
         for (i, v) in op.iter_nonzero() {
             self.act_on_basis(
                 result.copy(),
-                (coeff * v) % *p,
+                (coeff * v) % p,
                 op_degree,
                 i,
                 input_degree,
diff --git a/ext/crates/algebra/src/module/rpn.rs b/ext/crates/algebra/src/module/rpn.rs
index 2647a0afa..45677e350 100644
--- a/ext/crates/algebra/src/module/rpn.rs
+++ b/ext/crates/algebra/src/module/rpn.rs
@@ -184,7 +184,7 @@ where
     for<'a> &'a A: TryInto<&'a SteenrodAlgebra>,
 {
     pub fn new(algebra: Arc<A>, min: i32, max: Option<i32>, clear_bottom: bool) -> Self {
-        assert_eq!(*algebra.prime(), 2);
+        assert_eq!(algebra.prime(), 2);
         assert!(
             (&*algebra).try_into().is_ok(),
             "Real Projective Space only supports Steenrod Algebra"
diff --git a/ext/crates/algebra/src/module/tensor_module.rs b/ext/crates/algebra/src/module/tensor_module.rs
index 8c385b9f2..834bbb740 100644
--- a/ext/crates/algebra/src/module/tensor_module.rs
+++ b/ext/crates/algebra/src/module/tensor_module.rs
@@ -129,7 +129,7 @@ where
                         }
                         result.add_tensor(
                             self.offset(output_degree, left_deg + op_deg_l),
-                            minus_one_to_the_n(*self.prime(), op_deg_r * left_deg),
+                            minus_one_to_the_n(self.prime(), op_deg_r * left_deg),
                             left_result.as_slice(),
                             right_result.as_slice(),
                         );
diff --git a/ext/crates/algebra/src/steenrod_evaluator.rs b/ext/crates/algebra/src/steenrod_evaluator.rs
index 55e2cd1c7..95e659cb5 100644
--- a/ext/crates/algebra/src/steenrod_evaluator.rs
+++ b/ext/crates/algebra/src/steenrod_evaluator.rs
@@ -2,7 +2,7 @@ use crate::algebra::adem_algebra::AdemBasisElement;
 use crate::algebra::{AdemAlgebra, Algebra, MilnorAlgebra};
 use crate::milnor_algebra::{MilnorBasisElement, PPartEntry};
 use crate::steenrod_parser::*;
-use fp::prime::ValidPrime;
+use fp::prime::{Prime, ValidPrime};
 use fp::vector::FpVector;
 
 use anyhow::anyhow;
@@ -24,14 +24,14 @@ impl SteenrodEvaluator {
     pub fn milnor_to_adem(&self, result: &mut FpVector, coeff: u32, degree: i32, input: &FpVector) {
         let p = self.prime();
         for (i, v) in input.iter_nonzero() {
-            self.milnor_to_adem_on_basis(result, (coeff * v) % *p, degree, i);
+            self.milnor_to_adem_on_basis(result, (coeff * v) % p, degree, i);
         }
     }
 
     pub fn adem_to_milnor(&self, result: &mut FpVector, coeff: u32, degree: i32, input: &FpVector) {
         let p = self.prime();
         for (i, v) in input.iter_nonzero() {
-            self.adem_to_milnor_on_basis(result, (coeff * v) % *p, degree, i);
+            self.adem_to_milnor_on_basis(result, (coeff * v) % p, degree, i);
         }
     }
 
@@ -131,7 +131,7 @@ impl SteenrodEvaluator {
                     }
                 }
                 let mut result = FpVector::new(p, 1);
-                result.set_entry(0, x.rem_euclid(*p as i32) as u32);
+                result.set_entry(0, x.rem_euclid(p.as_i32()) as u32);
                 Ok((0, result))
             }
         }
@@ -274,7 +274,7 @@ impl SteenrodEvaluator {
             degree,
             bocksteins: 0,
             ps: t,
-            p_or_sq: *p != 2,
+            p_or_sq: p != 2,
         });
         let mut tmp_vector_a = FpVector::new(p, dim);
         self.adem_to_milnor_on_basis(&mut tmp_vector_a, 1, degree, t_idx);
@@ -315,19 +315,19 @@ impl SteenrodEvaluator {
             } else {
                 0
             };
-            t[i] = p_part + ((elt.q_part >> (i + 1)) & 1) + *p * t[i + 1];
+            t[i] = p_part + ((elt.q_part >> (i + 1)) & 1) + p * t[i + 1];
         }
         let t_idx = self.adem.basis_element_to_index(&AdemBasisElement {
             degree,
             bocksteins: elt.q_part,
             ps: t,
-            p_or_sq: *p != 2,
+            p_or_sq: p != 2,
         });
         let mut tmp_vector_a = FpVector::new(p, dim);
         self.adem_to_milnor_on_basis(&mut tmp_vector_a, 1, degree, t_idx);
         assert!(tmp_vector_a.entry(idx) == 1);
         tmp_vector_a.set_entry(idx, 0);
-        tmp_vector_a.scale(*p - 1);
+        tmp_vector_a.scale(p - 1);
         self.milnor_to_adem(result, coeff, degree, &tmp_vector_a);
         result.add_basis_element(t_idx, coeff);
     }
diff --git a/ext/crates/fp/benches/criterion.rs b/ext/crates/fp/benches/criterion.rs
index e9082ed30..440aed14c 100644
--- a/ext/crates/fp/benches/criterion.rs
+++ b/ext/crates/fp/benches/criterion.rs
@@ -33,7 +33,7 @@ fn row_reductions(c: &mut Criterion) {
 fn random_vector(p: ValidPrime, dimension: usize) -> Vec<u32> {
     let mut result = Vec::with_capacity(dimension);
     let mut rng = rand::thread_rng();
-    result.resize_with(dimension, || rng.gen::<u32>() % *p);
+    result.resize_with(dimension, || rng.gen::<u32>() % p);
     result
 }
 
diff --git a/ext/crates/fp/benches/iai.rs b/ext/crates/fp/benches/iai.rs
index 3ec9ba287..a7bd8f27b 100644
--- a/ext/crates/fp/benches/iai.rs
+++ b/ext/crates/fp/benches/iai.rs
@@ -13,7 +13,7 @@ fn random_matrix(p: ValidPrime, dimension: usize) -> Matrix {
 fn random_vector(p: ValidPrime, dimension: usize) -> Vec<u32> {
     let mut result = Vec::with_capacity(dimension);
     let mut rng = rand::thread_rng();
-    result.resize_with(dimension, || rng.gen::<u32>() % *p);
+    result.resize_with(dimension, || rng.gen::<u32>() % p);
     result
 }
 
diff --git a/ext/crates/fp/build.rs b/ext/crates/fp/build.rs
index 9d25e5a1b..62b3f28c4 100644
--- a/ext/crates/fp/build.rs
+++ b/ext/crates/fp/build.rs
@@ -20,10 +20,10 @@ fn main() -> Result<(), Error> {
 
     let mut writer = ConstWriter::for_build("constants")?.finish_dependencies();
 
-    writer.add_raw("/// The number of primes that are supported.");
+    writer.add_raw("/// The number of primes that will be prepared at compile-time. This includes");
+    writer.add_raw("/// precomputing inverses, binomial coefficients, etc.");
     writer.add_value("NUM_PRIMES", "usize", num_primes);
-    writer.add_raw("/// The `MAX_PRIME`th prime number. Constructing a `ValidPrime` using any number larger than");
-    writer.add_raw("/// this value will cause a panic.");
+    writer.add_raw("/// The `NUM_PRIMES`th prime number.");
     writer.add_value("MAX_PRIME", "usize", max_prime);
     // `NOT_A_PRIME` is never used if odd-primes is disabled.
     writer.add_raw("#[allow(dead_code)]");
diff --git a/ext/crates/fp/src/limb.rs b/ext/crates/fp/src/limb.rs
index 41e5a5e7c..3f5334b73 100644
--- a/ext/crates/fp/src/limb.rs
+++ b/ext/crates/fp/src/limb.rs
@@ -2,7 +2,7 @@ use std::ops::Range;
 
 pub(crate) use crate::constants::Limb;
 
-use crate::{constants::BITS_PER_LIMB, prime::ValidPrime};
+use crate::{constants::BITS_PER_LIMB, prime::Prime};
 
 /// A struct containing the information required to access a specific entry in an array of `Limb`s.
 #[derive(Copy, Clone)]
@@ -12,51 +12,38 @@ pub(crate) struct LimbBitIndexPair {
 }
 
 /// Return the number of bits an element of $\mathbb{F}_P$ occupies in a limb.
-pub(crate) const fn bit_length(p: ValidPrime) -> usize {
-    let p = p.value();
+pub(crate) fn bit_length<P: Prime>(p: P) -> usize {
+    let p = p.as_u32() as u64;
     match p {
         2 => 1,
-        _ => (32 - (p * (p - 1)).leading_zeros()) as usize,
-    }
-}
-
-/// Return the number of bits an element of $\mathbb{F}_P$ occupies in a limb.
-pub(crate) const fn bit_length_const<const P: u32>() -> usize {
-    match P {
-        2 => 1,
-        _ => (32 - (P * (P - 1)).leading_zeros()) as usize,
+        _ => (BITS_PER_LIMB as u32 - (p * (p - 1)).leading_zeros()) as usize,
     }
 }
 
 /// If `l` is a limb of elements of $\\mathbb{F}_p$, then `l & bitmask::<P>()` is the value of the
 /// first entry of `l`.
-pub(crate) const fn bitmask<const P: u32>() -> Limb {
-    (1 << bit_length_const::<P>()) - 1
+pub(crate) fn bitmask<P: Prime>(p: P) -> Limb {
+    (1 << bit_length(p)) - 1
 }
 
 // this function is never called if `odd-primes` is disabled
 #[allow(dead_code)]
 /// The number of elements of $\\mathbb{F}_p$ that fit in a single limb.
-pub(crate) const fn entries_per_limb(p: ValidPrime) -> usize {
+pub(crate) fn entries_per_limb<P: Prime>(p: P) -> usize {
     BITS_PER_LIMB / bit_length(p)
 }
 
-/// The number of elements of $\\mathbb{F}_p$ that fit in a single limb.
-pub(crate) const fn entries_per_limb_const<const P: u32>() -> usize {
-    BITS_PER_LIMB / bit_length_const::<P>()
-}
-
-pub(crate) const fn limb_bit_index_pair<const P: u32>(idx: usize) -> LimbBitIndexPair {
+pub(crate) fn limb_bit_index_pair<P: Prime>(p: P, idx: usize) -> LimbBitIndexPair {
     LimbBitIndexPair {
-        limb: idx / entries_per_limb_const::<P>(),
-        bit_index: (idx % entries_per_limb_const::<P>() * bit_length_const::<P>()),
+        limb: idx / entries_per_limb(p),
+        bit_index: (idx % entries_per_limb(p) * bit_length(p)),
     }
 }
 
 /// Return the `Limb` whose `i`th entry is `limb_a[i] + coeff * limb_b[i]` mod P. Both `limb_a` and
 /// `limb_b` are assumed to be reduced.
-pub(crate) const fn add<const P: u32>(limb_a: Limb, limb_b: Limb, coeff: u32) -> Limb {
-    if P == 2 {
+pub(crate) fn add<P: Prime>(p: P, limb_a: Limb, limb_b: Limb, coeff: u32) -> Limb {
+    if p == 2 {
         limb_a ^ (coeff as Limb * limb_b)
     } else {
         limb_a + (coeff as Limb) * limb_b
@@ -66,8 +53,8 @@ pub(crate) const fn add<const P: u32>(limb_a: Limb, limb_b: Limb, coeff: u32) ->
 /// Return the `Limb` whose entries are the entries of `limb` reduced modulo `P`.
 ///
 /// Contributed by Robert Burklund.
-pub(crate) fn reduce<const P: u32>(limb: Limb) -> Limb {
-    match P {
+pub(crate) fn reduce<P: Prime>(p: P, limb: Limb) -> Limb {
+    match p.as_u32() {
         2 => limb,
         3 => {
             // Set top bit to 1 in every limb
@@ -91,14 +78,14 @@ pub(crate) fn reduce<const P: u32>(limb: Limb) -> Limb {
             let d = m & BOTTOM_THREE_BITS;
             d + c - BOTTOM_TWO_BITS
         }
-        _ => pack::<_, P>(unpack::<P>(limb).map(|x| x % P)),
+        _ => pack(p, unpack(p, limb).map(|x| (x % (p.as_u32() as u64)) as u32)),
     }
 }
 
 /// Check whether or not a limb is reduced, i.e. whether every entry is a value in the range `0..P`.
 /// This is currently **not** faster than calling [`reduce`] directly.
-pub(crate) fn is_reduced<const P: u32>(limb: Limb) -> bool {
-    limb == reduce::<P>(limb)
+pub(crate) fn is_reduced<P: Prime>(p: P, limb: Limb) -> bool {
+    limb == reduce(p, limb)
 }
 
 /// Given an interator of `u32`'s, pack all of them into a single limb in order.
@@ -107,8 +94,8 @@ pub(crate) fn is_reduced<const P: u32>(limb: Limb) -> bool {
 ///  - The values of the iterator fit into a single limb
 ///
 /// If these assumptions are violated, the result will be nonsense.
-pub(crate) fn pack<T: Iterator<Item = u32>, const P: u32>(entries: T) -> Limb {
-    let bit_length = bit_length_const::<P>();
+pub(crate) fn pack<T: Iterator<Item = u32>, P: Prime>(p: P, entries: T) -> Limb {
+    let bit_length = bit_length(p);
     let mut result: Limb = 0;
     let mut shift = 0;
     for entry in entries {
@@ -118,34 +105,36 @@ pub(crate) fn pack<T: Iterator<Item = u32>, const P: u32>(entries: T) -> Limb {
     result
 }
 
-/// Give an iterator over the entries of `limb`.
-pub(crate) fn unpack<const P: u32>(mut limb: Limb) -> impl Iterator<Item = u32> {
-    let entries = entries_per_limb_const::<P>();
-    let bit_length = bit_length_const::<P>();
-    let bit_mask = bitmask::<P>();
+/// Give an iterator over the entries of `limb`. We return `u64` instead of `u32` because the
+/// entries might not be reduced. In general we can only guarantee that the entries are less than p
+/// * (p - 1), so they might not fit into a single `u32`.
+pub(crate) fn unpack<P: Prime>(p: P, mut limb: Limb) -> impl Iterator<Item = u64> {
+    let entries = entries_per_limb(p);
+    let bit_length = bit_length(p);
+    let bit_mask = bitmask(p);
 
     (0..entries).map(move |_| {
-        let result = (limb & bit_mask) as u32;
+        let result = limb & bit_mask;
         limb >>= bit_length;
         result
     })
 }
 
 /// Return the number of limbs required to hold `dim` entries.
-pub(crate) const fn number<const P: u32>(dim: usize) -> usize {
+pub(crate) fn number<P: Prime>(p: P, dim: usize) -> usize {
     if dim == 0 {
         0
     } else {
-        limb_bit_index_pair::<P>(dim - 1).limb + 1
+        limb_bit_index_pair(p, dim - 1).limb + 1
     }
 }
 
 /// Return the `Range<usize>` starting at the index of the limb containing the `start`th entry, and
 /// ending at the index of the limb containing the `end`th entry (including the latter).
-pub(crate) const fn range<const P: u32>(start: usize, end: usize) -> Range<usize> {
-    let min = limb_bit_index_pair::<P>(start).limb;
+pub(crate) fn range<P: Prime>(p: P, start: usize, end: usize) -> Range<usize> {
+    let min = limb_bit_index_pair(p, start).limb;
     let max = if end > 0 {
-        limb_bit_index_pair::<P>(end - 1).limb + 1
+        limb_bit_index_pair(p, end - 1).limb + 1
     } else {
         0
     };
@@ -169,8 +158,8 @@ pub(crate) fn sign_rule(mut target: Limb, mut source: Limb) -> u32 {
 }
 
 /// Return either `Some(sum)` if no carries happen in the limb, or `None` if some carry does happen.
-pub(crate) fn truncate<const P: u32>(sum: Limb) -> Option<Limb> {
-    if is_reduced::<P>(sum) {
+pub(crate) fn truncate<P: Prime>(p: P, sum: Limb) -> Option<Limb> {
+    if is_reduced(p, sum) {
         Some(sum)
     } else {
         None
diff --git a/ext/crates/fp/src/matrix/m4ri.rs b/ext/crates/fp/src/matrix/m4ri.rs
index 8837823a4..6a61ae0d7 100644
--- a/ext/crates/fp/src/matrix/m4ri.rs
+++ b/ext/crates/fp/src/matrix/m4ri.rs
@@ -1,6 +1,7 @@
 use crate::constants::BITS_PER_LIMB;
 use crate::limb::{number, Limb};
 use crate::matrix::Matrix;
+use crate::prime::P2;
 use crate::simd;
 
 use itertools::Itertools;
@@ -35,7 +36,7 @@ pub(crate) struct M4riTable {
 impl M4riTable {
     /// Create a table with space for `k` vectors, each with `cols` columns.
     pub fn new(k: usize, cols: usize) -> Self {
-        let num_limbs = number::<2>(cols);
+        let num_limbs = number(P2, cols);
         Self {
             rows: Vec::with_capacity(k),
             columns: Vec::with_capacity(k),
diff --git a/ext/crates/fp/src/matrix/matrix_inner.rs b/ext/crates/fp/src/matrix/matrix_inner.rs
index 7a703bdba..b4b46d160 100644
--- a/ext/crates/fp/src/matrix/matrix_inner.rs
+++ b/ext/crates/fp/src/matrix/matrix_inner.rs
@@ -393,7 +393,7 @@ impl Matrix {
             return;
         }
         let (target, source) = self.split_borrow(target, source);
-        target.add_offset(source, *prime - coef, pivot_column);
+        target.add_offset(source, prime - coef, pivot_column);
     }
 
     /// A version of [`Matrix::row_op`] without the zero assumption.
@@ -413,7 +413,7 @@ impl Matrix {
             return;
         }
         let (target, source) = self.split_borrow(target, source);
-        target.add(source, *prime - coef);
+        target.add(source, prime - coef);
     }
 
     /// Mutably borrows `x[i]` and `x[j]`.
@@ -516,7 +516,7 @@ impl Matrix {
 
         let mut empty_rows = Vec::with_capacity(self.rows());
 
-        if *self.p == 2 {
+        if self.p == 2 {
             // the m4ri C library uses a similar formula but with a hard cap of 7 instead of 8
             let k = std::cmp::min(8, crate::prime::log2(1 + self.rows()) * 3 / 4);
             let mut table = M4riTable::new(k, self.columns());
@@ -882,7 +882,7 @@ impl Matrix {
         for i in 0..input.len() {
             result.add(
                 self.vectors[i].as_slice(),
-                (coeff * input.entry(i)) % *self.p,
+                (coeff * input.entry(i)) % self.p,
             );
         }
     }
@@ -919,7 +919,7 @@ impl std::ops::Mul for &Matrix {
 impl std::ops::MulAssign<u32> for Matrix {
     fn mul_assign(&mut self, rhs: u32) {
         #[allow(clippy::suspicious_op_assign_impl)]
-        let rhs = rhs % *self.p;
+        let rhs = rhs % self.p;
         for row in self.iter_mut() {
             row.scale(rhs);
         }
diff --git a/ext/crates/fp/src/matrix/quasi_inverse.rs b/ext/crates/fp/src/matrix/quasi_inverse.rs
index 276015114..fb8d5c012 100644
--- a/ext/crates/fp/src/matrix/quasi_inverse.rs
+++ b/ext/crates/fp/src/matrix/quasi_inverse.rs
@@ -138,7 +138,7 @@ impl QuasiInverse {
                 }
             }
             if c != 0 {
-                target.add(self.preimage[row].as_slice(), (coeff * c) % *p);
+                target.add(self.preimage[row].as_slice(), (coeff * c) % p);
             }
             row += 1;
         }
diff --git a/ext/crates/fp/src/matrix/subquotient.rs b/ext/crates/fp/src/matrix/subquotient.rs
index 0aeca4c32..c1302239b 100644
--- a/ext/crates/fp/src/matrix/subquotient.rs
+++ b/ext/crates/fp/src/matrix/subquotient.rs
@@ -52,7 +52,7 @@ impl Subquotient {
             if c != 0 {
                 elt.add(
                     self.gens.matrix.row(self.gens.pivots()[i] as usize),
-                    ((*elt.prime() - 1) * c) % *elt.prime(),
+                    ((elt.prime() - 1) * c) % elt.prime(),
                 );
             }
         }
diff --git a/ext/crates/fp/src/matrix/subspace.rs b/ext/crates/fp/src/matrix/subspace.rs
index 488a684b6..80997da74 100644
--- a/ext/crates/fp/src/matrix/subspace.rs
+++ b/ext/crates/fp/src/matrix/subspace.rs
@@ -129,7 +129,7 @@ impl Subspace {
         for (col, row) in iter {
             let c = vector.as_slice().entry(col);
             if c != 0 {
-                vector.add(row, *p - c);
+                vector.add(row, p - c);
             }
         }
     }
diff --git a/ext/crates/fp/src/prime.rs b/ext/crates/fp/src/prime.rs
deleted file mode 100644
index a5ff1e6b2..000000000
--- a/ext/crates/fp/src/prime.rs
+++ /dev/null
@@ -1,593 +0,0 @@
-use serde::{de::Error, Deserialize, Deserializer, Serialize, Serializer};
-use std::convert::TryFrom;
-
-use crate::constants::{
-    BINOMIAL4_TABLE, BINOMIAL4_TABLE_SIZE, BINOMIAL_TABLE, INVERSE_TABLE, PRIME_TO_INDEX_MAP,
-};
-#[allow(unused_imports)]
-use crate::constants::{MAX_PRIME, NOT_A_PRIME, NUM_PRIMES};
-
-pub const TWO: ValidPrime = ValidPrime::new(2);
-
-#[derive(Copy, Clone, Debug, PartialEq, Eq)]
-pub struct ValidPrime {
-    #[cfg(feature = "odd-primes")]
-    p: u32,
-}
-
-impl ValidPrime {
-    pub const fn new(p: u32) -> Self {
-        assert!(is_valid_prime(p));
-
-        #[cfg(feature = "odd-primes")]
-        {
-            Self { p }
-        }
-
-        #[cfg(not(feature = "odd-primes"))]
-        {
-            Self {}
-        }
-    }
-
-    /// Get the underlying prime. This is the same function as deref but
-    /// 1. This is a const fn
-    /// 2. This does not inform the compiler about properties of p via unreachable_unchecked.
-    ///
-    /// Use this function in a const context where you would expect it to be evaluated at
-    /// compile-time.
-    pub const fn value(&self) -> u32 {
-        #[cfg(feature = "odd-primes")]
-        {
-            self.p
-        }
-        #[cfg(not(feature = "odd-primes"))]
-        {
-            2
-        }
-    }
-}
-
-impl std::ops::Deref for ValidPrime {
-    type Target = u32;
-
-    #[cfg(not(feature = "odd-primes"))]
-    fn deref(&self) -> &Self::Target {
-        &2
-    }
-
-    #[cfg(feature = "odd-primes")]
-    fn deref(&self) -> &Self::Target {
-        let p = self.p;
-        unsafe {
-            if !is_valid_prime(p) || p == 0 || PRIME_TO_INDEX_MAP[p as usize] >= NUM_PRIMES {
-                std::hint::unreachable_unchecked()
-            }
-        }
-        &self.p
-    }
-}
-
-#[derive(Debug)]
-pub struct InvalidPrimeError(u32);
-
-impl std::fmt::Display for InvalidPrimeError {
-    fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result {
-        write!(f, "Invalid prime: {}", self.0)
-    }
-}
-
-impl std::error::Error for InvalidPrimeError {}
-
-impl TryFrom<u32> for ValidPrime {
-    type Error = InvalidPrimeError;
-
-    fn try_from(p: u32) -> Result<Self, InvalidPrimeError> {
-        if is_valid_prime(p) {
-            Ok(Self::new(p))
-        } else {
-            Err(InvalidPrimeError(p))
-        }
-    }
-}
-
-impl std::str::FromStr for ValidPrime {
-    type Err = String;
-
-    fn from_str(s: &str) -> Result<Self, String> {
-        let p: u32 = s.parse::<u32>().map_err(|s| s.to_string())?;
-        Self::try_from(p).map_err(|s| s.to_string())
-    }
-}
-
-impl From<ValidPrime> for u32 {
-    fn from(p: ValidPrime) -> u32 {
-        *p
-    }
-}
-
-impl std::fmt::Display for ValidPrime {
-    fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result {
-        (**self).fmt(f)
-    }
-}
-
-impl Serialize for ValidPrime {
-    fn serialize<S>(&self, serializer: S) -> Result<S::Ok, S::Error>
-    where
-        S: Serializer,
-    {
-        (**self).serialize(serializer)
-    }
-}
-
-impl<'de> Deserialize<'de> for ValidPrime {
-    fn deserialize<D>(deserializer: D) -> Result<Self, D::Error>
-    where
-        D: Deserializer<'de>,
-    {
-        let p: u32 = u32::deserialize(deserializer)?;
-        ValidPrime::try_from(p).map_err(D::Error::custom)
-    }
-}
-
-#[cfg(feature = "odd-primes")]
-pub const fn is_valid_prime(p: u32) -> bool {
-    (p as usize) <= MAX_PRIME && PRIME_TO_INDEX_MAP[p as usize] != NOT_A_PRIME
-}
-
-#[cfg(not(feature = "odd-primes"))]
-pub const fn is_valid_prime(p: u32) -> bool {
-    p == 2
-}
-
-// Uses a the lookup table we initialized.
-pub fn inverse(p: ValidPrime, k: u32) -> u32 {
-    assert!(k > 0 && k < *p);
-    // LLVM doesn't understand the inequality is transitive
-    unsafe { *INVERSE_TABLE[PRIME_TO_INDEX_MAP[*p as usize]].get_unchecked(k as usize) }
-}
-
-pub const fn minus_one_to_the_n(p: u32, i: i32) -> u32 {
-    if i % 2 == 0 {
-        1
-    } else {
-        p - 1
-    }
-}
-
-/// This uses a lookup table for n choose k when n and k are both less than p.
-/// Lucas's theorem reduces general binomial coefficients to this case.
-///
-/// Calling this function safely requires that `k, n < p`.  These invariants are often known
-/// apriori because k and n are obtained by reducing mod p, so it is better to expose an unsafe
-/// interface that avoids these checks.
-unsafe fn direct_binomial(p: ValidPrime, n: usize, k: usize) -> u32 {
-    *BINOMIAL_TABLE
-        .get_unchecked(PRIME_TO_INDEX_MAP[*p as usize])
-        .get_unchecked(n)
-        .get_unchecked(k)
-}
-
-/// Computes b^e.
-pub const fn integer_power(mut b: u32, mut e: u32) -> u32 {
-    let mut result: u32 = 1;
-    while e > 0 {
-        // b is b^{2^i}
-        // if the current bit of e is odd, mutliply b^{2^i} into result.
-        if e & 1 == 1 {
-            result *= b;
-        }
-        b *= b;
-        e >>= 1;
-    }
-    result
-}
-
-/// Compute b^e mod p.
-pub const fn power_mod(p: u32, mut b: u32, mut e: u32) -> u32 {
-    assert!(p > 0);
-    let mut result: u32 = 1;
-    while e > 0 {
-        if (e & 1) == 1 {
-            result = (result * b) % p;
-        }
-        b = (b * b) % p;
-        e >>= 1;
-    }
-    result
-}
-
-/// Compute the base 2 log of a number, rounded down to the nearest integer.
-///
-/// # Example
-/// ```
-/// # use fp::prime::log2;
-/// assert_eq!(0, log2(0b1));
-/// assert_eq!(1, log2(0b10));
-/// assert_eq!(1, log2(0b11));
-/// assert_eq!(3, log2(0b1011));
-/// ```
-pub const fn log2(n: usize) -> usize {
-    std::mem::size_of::<usize>() * 8 - 1 - n.leading_zeros() as usize
-}
-
-/// Discrete log base p of n.
-pub const fn logp(p: u32, mut n: u32) -> u32 {
-    let mut result = 0u32;
-    while n > 0 {
-        n /= p;
-        result += 1;
-    }
-    result
-}
-
-/// Factor $n$ as $p^k m$. Returns $(k, m)$.
-pub const fn factor_pk(p: u32, mut n: u32) -> (u32, u32) {
-    if n == 0 {
-        return (0, 0);
-    }
-    let mut k = 0;
-    while n % p == 0 {
-        n /= p;
-        k += 1;
-    }
-    (k, n)
-}
-/// A number satisfying the Binomial trait supports computing various binomial coefficients. This
-/// is implemented using a macro, since the implementation for all types is syntactically the same.
-pub trait Binomial: Sized {
-    /// mod 2 multinomial coefficient
-    fn multinomial2(k: &[Self]) -> Self;
-
-    /// mod 2 binomial coefficient n choose k
-    fn binomial2(n: Self, k: Self) -> Self;
-
-    /// Binomial coefficients mod 4. We pre-compute the coefficients for small values of n. For large
-    /// n, we recursively use the fact that if n = 2^k + l, l < 2^k, then
-    ///
-    ///    n choose r = l choose r + 2 (l choose (r - 2^{k - 1})) + (l choose (r - 2^k))
-    ///
-    /// This is easy to verify using the fact that
-    ///
-    ///    (x + y)^{2^k} = x^{2^k} + 2 x^{2^{k - 1}} y^{2^{k - 1}} + y^{2^k}
-    ///
-    fn binomial4(n: Self, k: Self) -> Self;
-
-    /// Compute binomial coefficients mod 4 using the recursion relation in the documentation of
-    /// [Binomial::binomial4]. This calls into binomial4 instead of binomial4_rec. The main purpose
-    /// of this is to separate out the logic for testing.
-    fn binomial4_rec(n: Self, k: Self) -> Self;
-
-    /// Computes the multinomial coefficient mod p using Lucas' theorem. This modifies the
-    /// underlying list. For p = 2 it is more efficient to use multinomial2
-    fn multinomial_odd(p: ValidPrime, l: &mut [Self]) -> Self;
-
-    /// Compute odd binomial coefficients mod p, where p is odd. For p = 2 it is more efficient to
-    /// use binomial2
-    fn binomial_odd(p: ValidPrime, n: Self, k: Self) -> Self;
-
-    /// Checks whether n choose k is zero mod p. Since we don't have to compute the value, this is
-    /// faster than binomial_odd.
-    fn binomial_odd_is_zero(p: ValidPrime, n: Self, k: Self) -> bool;
-
-    /// Multinomial coefficient of the list l
-    fn multinomial(p: ValidPrime, l: &mut [Self]) -> Self {
-        if *p == 2 {
-            Self::multinomial2(l)
-        } else {
-            Self::multinomial_odd(p, l)
-        }
-    }
-
-    /// Binomial coefficient n choose k.
-    fn binomial(p: ValidPrime, n: Self, k: Self) -> Self {
-        if *p == 2 {
-            Self::binomial2(n, k)
-        } else {
-            Self::binomial_odd(p, n, k)
-        }
-    }
-}
-
-macro_rules! impl_binomial {
-    ($T:ty) => {
-        impl Binomial for $T {
-            #[inline]
-            fn multinomial2(l: &[Self]) -> Self {
-                let mut bit_or: Self = 0;
-                let mut sum: Self = 0;
-                for &e in l {
-                    sum += e;
-                    bit_or |= e;
-                }
-                if bit_or == sum {
-                    1
-                } else {
-                    0
-                }
-            }
-
-            #[inline]
-            fn binomial2(n: Self, k: Self) -> Self {
-                if n < k {
-                    0
-                } else if (n - k) & k == 0 {
-                    1
-                } else {
-                    0
-                }
-            }
-            #[inline]
-            fn multinomial_odd(p_: ValidPrime, l: &mut [Self]) -> Self {
-                let p = *p_ as Self;
-
-                let mut n: Self = l.iter().sum();
-                if n == 0 {
-                    return 1;
-                }
-                let mut answer = 1;
-
-                while n > 0 {
-                    let mut multi: Self = 1;
-
-                    let total_entry = n % p;
-                    n /= p;
-
-                    let mut partial_sum: Self = l[0] % p;
-                    l[0] /= p;
-
-                    for ll in l.iter_mut().skip(1) {
-                        let entry = *ll % p;
-                        *ll /= p;
-
-                        partial_sum += entry;
-                        if partial_sum > total_entry {
-                            // This early return is necessary because direct_binomial only works when
-                            // partial_sum < 19
-                            return 0;
-                        }
-                        // This is safe because p < 20, partial_sum <= total_entry < p and entry < p.
-                        multi *=
-                            unsafe { direct_binomial(p_, partial_sum as usize, entry as usize) }
-                                as Self;
-                        multi %= p;
-                    }
-                    answer *= multi;
-                    answer %= p;
-                }
-                answer
-            }
-
-            #[inline]
-            fn binomial_odd(p_: ValidPrime, mut n: Self, mut k: Self) -> Self {
-                let p = *p_ as Self;
-
-                // We have both signed and unsigned types
-                #[allow(unused_comparisons)]
-                if n < k || k < 0 {
-                    return 0;
-                }
-
-                let mut answer = 1;
-
-                while n > 0 {
-                    // This is safe because p < 20 and anything mod p is < p.
-                    answer *=
-                        unsafe { direct_binomial(p_, (n % p) as usize, (k % p) as usize) } as Self;
-                    answer %= p;
-                    n /= p;
-                    k /= p;
-                }
-                answer
-            }
-
-            #[inline]
-            fn binomial_odd_is_zero(p: ValidPrime, mut n: Self, mut k: Self) -> bool {
-                let p = *p as Self;
-
-                while n > 0 {
-                    if n % p < k % p {
-                        return true;
-                    }
-                    n /= p;
-                    k /= p;
-                }
-                false
-            }
-            fn binomial4(n: Self, j: Self) -> Self {
-                if (n as usize) < BINOMIAL4_TABLE_SIZE {
-                    return BINOMIAL4_TABLE[n as usize][j as usize] as Self;
-                }
-                if (n - j) & j == 0 {
-                    // Answer is odd
-                    Self::binomial4_rec(n, j)
-                } else if (n - j).count_ones() + j.count_ones() - n.count_ones() == 1 {
-                    2
-                } else {
-                    0
-                }
-            }
-
-            #[inline]
-            fn binomial4_rec(n: Self, j: Self) -> Self {
-                let k = (std::mem::size_of::<Self>() * 8) as u32 - n.leading_zeros() - 1;
-                let l = n - (1 << k);
-                let mut ans = 0;
-                if j <= l {
-                    ans += Self::binomial4(l, j)
-                }
-                let pow = 1 << (k - 1);
-                if pow <= j && j <= l + pow {
-                    ans += 2 * Self::binomial2(l, j - pow);
-                }
-                if j >= (1 << k) {
-                    ans += Self::binomial4(l, j - (1 << k));
-                }
-                ans % 4
-            }
-        }
-    };
-}
-
-impl_binomial!(u32);
-impl_binomial!(u16);
-impl_binomial!(i32);
-
-pub struct BitflagIterator {
-    remaining: u8,
-    flag: u64,
-}
-
-impl BitflagIterator {
-    pub fn new(flag: u64) -> Self {
-        Self {
-            remaining: u8::max_value(),
-            flag,
-        }
-    }
-
-    pub fn new_fixed_length(flag: u64, remaining: usize) -> Self {
-        assert!(remaining <= 64);
-        let remaining = remaining as u8;
-        Self { remaining, flag }
-    }
-
-    pub fn set_bit_iterator(flag: u64) -> impl Iterator<Item = usize> {
-        Self::new(flag)
-            .enumerate()
-            .filter_map(|(idx, v)| if v { Some(idx) } else { None })
-    }
-}
-
-impl Iterator for BitflagIterator {
-    type Item = bool;
-    fn next(&mut self) -> Option<Self::Item> {
-        if self.remaining > 64 && self.flag == 0 || self.remaining == 0 {
-            None
-        } else {
-            self.remaining -= 1;
-            let result = self.flag & 1 == 1;
-            self.flag >>= 1;
-            Some(result)
-        }
-    }
-}
-
-/// Iterates through all numbers with the same number of bits. It may panic or return nonsense
-/// after all valid values are returned.
-pub struct BinomialIterator {
-    value: u32,
-}
-
-impl BinomialIterator {
-    pub fn new(n: usize) -> Self {
-        Self {
-            value: (1 << n) - 1,
-        }
-    }
-}
-
-impl Iterator for BinomialIterator {
-    type Item = u32;
-    fn next(&mut self) -> Option<Self::Item> {
-        let v = self.value;
-        let c = v & v.wrapping_neg();
-        let r = v + c;
-        let n = (r ^ v).wrapping_shr(2 + v.trailing_zeros());
-        self.value = r | n;
-        Some(v)
-    }
-}
-
-#[cfg(test)]
-mod tests {
-    use super::*;
-    use crate::constants::PRIMES;
-
-    #[test]
-    fn inverse_test() {
-        for &p in PRIMES.iter() {
-            let p = ValidPrime::new(p);
-            for k in 1..*p {
-                assert_eq!((inverse(p, k) * k) % *p, 1);
-            }
-        }
-    }
-
-    #[test]
-    fn binomial_test() {
-        let entries = [[2, 2, 1, 0], [2, 3, 1, 1], [3, 1090, 730, 1], [7, 3, 2, 3]];
-
-        for entry in &entries {
-            assert_eq!(
-                entry[3],
-                u32::binomial(ValidPrime::new(entry[0]), entry[1], entry[2])
-            );
-        }
-    }
-
-    #[test]
-    fn binomial_vs_monomial() {
-        for &p in &[2, 3, 5, 7, 11] {
-            let p = ValidPrime::new(p);
-            for l in 0..20 {
-                for m in 0..20 {
-                    assert_eq!(u32::binomial(p, l + m, m), u32::multinomial(p, &mut [l, m]))
-                }
-            }
-        }
-    }
-
-    fn binomial_full(n: u32, j: u32) -> u32 {
-        let mut res = 1;
-        for k in j + 1..=n {
-            res *= k;
-        }
-        for k in 1..=(n - j) {
-            res /= k;
-        }
-        res
-    }
-
-    #[test]
-    fn binomial_cmp() {
-        for n in 0..12 {
-            for j in 0..=n {
-                let ans = binomial_full(n, j);
-                for &p in &[2, 3, 5, 7, 11] {
-                    assert_eq!(
-                        u32::binomial(ValidPrime::new(p), n, j),
-                        ans % p,
-                        "{n} choose {j} mod {p}"
-                    );
-                }
-                assert_eq!(u32::binomial4(n, j), ans % 4, "{n} choose {j} mod 4");
-                // binomial4_rec is only called on large n. It does not handle the n = 0, 1 cases
-                // correctly.
-                if n > 1 {
-                    assert_eq!(
-                        u32::binomial4_rec(n, j),
-                        ans % 4,
-                        "{n} choose {j} mod 4 rec"
-                    );
-                }
-            }
-        }
-    }
-
-    #[test]
-    fn binomial_iterator() {
-        let mut iter = BinomialIterator::new(4);
-        assert_eq!(iter.next(), Some(0b1111));
-        assert_eq!(iter.next(), Some(0b10111));
-        assert_eq!(iter.next(), Some(0b11011));
-        assert_eq!(iter.next(), Some(0b11101));
-        assert_eq!(iter.next(), Some(0b11110));
-        assert_eq!(iter.next(), Some(0b100111));
-        assert_eq!(iter.next(), Some(0b101011));
-        assert_eq!(iter.next(), Some(0b101101));
-        assert_eq!(iter.next(), Some(0b101110));
-        assert_eq!(iter.next(), Some(0b110011));
-    }
-}
diff --git a/ext/crates/fp/src/prime/binomial.rs b/ext/crates/fp/src/prime/binomial.rs
new file mode 100644
index 000000000..4fcf3fe6d
--- /dev/null
+++ b/ext/crates/fp/src/prime/binomial.rs
@@ -0,0 +1,219 @@
+use crate::{
+    constants::{BINOMIAL4_TABLE, BINOMIAL4_TABLE_SIZE, BINOMIAL_TABLE},
+    PRIME_TO_INDEX_MAP,
+};
+
+use super::{Prime, ValidPrime};
+
+/// This uses a lookup table for n choose k when n and k are both less than p.
+/// Lucas's theorem reduces general binomial coefficients to this case.
+///
+/// Calling this function safely requires that `k, n < p`.  These invariants are often known
+/// apriori because k and n are obtained by reducing mod p, so it is better to expose an unsafe
+/// interface that avoids these checks.
+unsafe fn direct_binomial(p: ValidPrime, n: usize, k: usize) -> u32 {
+    *BINOMIAL_TABLE
+        .get_unchecked(PRIME_TO_INDEX_MAP[p.as_usize()])
+        .get_unchecked(n)
+        .get_unchecked(k)
+}
+
+/// A number satisfying the Binomial trait supports computing various binomial coefficients. This
+/// is implemented using a macro, since the implementation for all types is syntactically the same.
+pub trait Binomial: Sized {
+    /// mod 2 multinomial coefficient
+    fn multinomial2(k: &[Self]) -> Self;
+
+    /// mod 2 binomial coefficient n choose k
+    fn binomial2(n: Self, k: Self) -> Self;
+
+    /// Binomial coefficients mod 4. We pre-compute the coefficients for small values of n. For large
+    /// n, we recursively use the fact that if n = 2^k + l, l < 2^k, then
+    ///
+    ///    n choose r = l choose r + 2 (l choose (r - 2^{k - 1})) + (l choose (r - 2^k))
+    ///
+    /// This is easy to verify using the fact that
+    ///
+    ///    (x + y)^{2^k} = x^{2^k} + 2 x^{2^{k - 1}} y^{2^{k - 1}} + y^{2^k}
+    ///
+    fn binomial4(n: Self, k: Self) -> Self;
+
+    /// Compute binomial coefficients mod 4 using the recursion relation in the documentation of
+    /// [Binomial::binomial4]. This calls into binomial4 instead of binomial4_rec. The main purpose
+    /// of this is to separate out the logic for testing.
+    fn binomial4_rec(n: Self, k: Self) -> Self;
+
+    /// Computes the multinomial coefficient mod p using Lucas' theorem. This modifies the
+    /// underlying list. For p = 2 it is more efficient to use multinomial2
+    fn multinomial_odd(p: ValidPrime, l: &mut [Self]) -> Self;
+
+    /// Compute odd binomial coefficients mod p, where p is odd. For p = 2 it is more efficient to
+    /// use binomial2
+    fn binomial_odd(p: ValidPrime, n: Self, k: Self) -> Self;
+
+    /// Checks whether n choose k is zero mod p. Since we don't have to compute the value, this is
+    /// faster than binomial_odd.
+    fn binomial_odd_is_zero(p: ValidPrime, n: Self, k: Self) -> bool;
+
+    /// Multinomial coefficient of the list l
+    fn multinomial(p: ValidPrime, l: &mut [Self]) -> Self {
+        if p == 2 {
+            Self::multinomial2(l)
+        } else {
+            Self::multinomial_odd(p, l)
+        }
+    }
+
+    /// Binomial coefficient n choose k.
+    fn binomial(p: ValidPrime, n: Self, k: Self) -> Self {
+        if p == 2 {
+            Self::binomial2(n, k)
+        } else {
+            Self::binomial_odd(p, n, k)
+        }
+    }
+}
+
+macro_rules! impl_binomial {
+    ($T:ty) => {
+        impl Binomial for $T {
+            #[inline]
+            fn multinomial2(l: &[Self]) -> Self {
+                let mut bit_or: Self = 0;
+                let mut sum: Self = 0;
+                for &e in l {
+                    sum += e;
+                    bit_or |= e;
+                }
+                if bit_or == sum {
+                    1
+                } else {
+                    0
+                }
+            }
+
+            #[inline]
+            fn binomial2(n: Self, k: Self) -> Self {
+                if n < k {
+                    0
+                } else if (n - k) & k == 0 {
+                    1
+                } else {
+                    0
+                }
+            }
+            #[inline]
+            fn multinomial_odd(p_: ValidPrime, l: &mut [Self]) -> Self {
+                let p = p_.as_u32() as Self;
+
+                let mut n: Self = l.iter().sum();
+                if n == 0 {
+                    return 1;
+                }
+                let mut answer = 1;
+
+                while n > 0 {
+                    let mut multi: Self = 1;
+
+                    let total_entry = n % p;
+                    n /= p;
+
+                    let mut partial_sum: Self = l[0] % p;
+                    l[0] /= p;
+
+                    for ll in l.iter_mut().skip(1) {
+                        let entry = *ll % p;
+                        *ll /= p;
+
+                        partial_sum += entry;
+                        if partial_sum > total_entry {
+                            // This early return is necessary because direct_binomial only works when
+                            // partial_sum < 19
+                            return 0;
+                        }
+                        // This is safe because p < 20, partial_sum <= total_entry < p and entry < p.
+                        multi *=
+                            unsafe { direct_binomial(p_, partial_sum as usize, entry as usize) }
+                                as Self;
+                        multi %= p;
+                    }
+                    answer *= multi;
+                    answer %= p;
+                }
+                answer
+            }
+
+            #[inline]
+            fn binomial_odd(p_: ValidPrime, mut n: Self, mut k: Self) -> Self {
+                let p = p_.as_u32() as Self;
+
+                // We have both signed and unsigned types
+                #[allow(unused_comparisons)]
+                if n < k || k < 0 {
+                    return 0;
+                }
+
+                let mut answer = 1;
+
+                while n > 0 {
+                    // This is safe because p < 20 and anything mod p is < p.
+                    answer *=
+                        unsafe { direct_binomial(p_, (n % p) as usize, (k % p) as usize) } as Self;
+                    answer %= p;
+                    n /= p;
+                    k /= p;
+                }
+                answer
+            }
+
+            #[inline]
+            fn binomial_odd_is_zero(p: ValidPrime, mut n: Self, mut k: Self) -> bool {
+                let p = p.as_u32() as Self;
+
+                while n > 0 {
+                    if n % p < k % p {
+                        return true;
+                    }
+                    n /= p;
+                    k /= p;
+                }
+                false
+            }
+            fn binomial4(n: Self, j: Self) -> Self {
+                if (n as usize) < BINOMIAL4_TABLE_SIZE {
+                    return BINOMIAL4_TABLE[n as usize][j as usize] as Self;
+                }
+                if (n - j) & j == 0 {
+                    // Answer is odd
+                    Self::binomial4_rec(n, j)
+                } else if (n - j).count_ones() + j.count_ones() - n.count_ones() == 1 {
+                    2
+                } else {
+                    0
+                }
+            }
+
+            #[inline]
+            fn binomial4_rec(n: Self, j: Self) -> Self {
+                let k = (std::mem::size_of::<Self>() * 8) as u32 - n.leading_zeros() - 1;
+                let l = n - (1 << k);
+                let mut ans = 0;
+                if j <= l {
+                    ans += Self::binomial4(l, j)
+                }
+                let pow = 1 << (k - 1);
+                if pow <= j && j <= l + pow {
+                    ans += 2 * Self::binomial2(l, j - pow);
+                }
+                if j >= (1 << k) {
+                    ans += Self::binomial4(l, j - (1 << k));
+                }
+                ans % 4
+            }
+        }
+    };
+}
+
+impl_binomial!(u32);
+impl_binomial!(u16);
+impl_binomial!(i32);
diff --git a/ext/crates/fp/src/prime/iter.rs b/ext/crates/fp/src/prime/iter.rs
new file mode 100644
index 000000000..a8fe72654
--- /dev/null
+++ b/ext/crates/fp/src/prime/iter.rs
@@ -0,0 +1,65 @@
+pub struct BitflagIterator {
+    remaining: u8,
+    flag: u64,
+}
+
+impl BitflagIterator {
+    pub fn new(flag: u64) -> Self {
+        Self {
+            remaining: u8::max_value(),
+            flag,
+        }
+    }
+
+    pub fn new_fixed_length(flag: u64, remaining: usize) -> Self {
+        assert!(remaining <= 64);
+        let remaining = remaining as u8;
+        Self { remaining, flag }
+    }
+
+    pub fn set_bit_iterator(flag: u64) -> impl Iterator<Item = usize> {
+        Self::new(flag)
+            .enumerate()
+            .filter_map(|(idx, v)| if v { Some(idx) } else { None })
+    }
+}
+
+impl Iterator for BitflagIterator {
+    type Item = bool;
+    fn next(&mut self) -> Option<Self::Item> {
+        if self.remaining > 64 && self.flag == 0 || self.remaining == 0 {
+            None
+        } else {
+            self.remaining -= 1;
+            let result = self.flag & 1 == 1;
+            self.flag >>= 1;
+            Some(result)
+        }
+    }
+}
+
+/// Iterates through all numbers with the same number of bits. It may panic or return nonsense
+/// after all valid values are returned.
+pub struct BinomialIterator {
+    value: u32,
+}
+
+impl BinomialIterator {
+    pub fn new(n: usize) -> Self {
+        Self {
+            value: (1 << n) - 1,
+        }
+    }
+}
+
+impl Iterator for BinomialIterator {
+    type Item = u32;
+    fn next(&mut self) -> Option<Self::Item> {
+        let v = self.value;
+        let c = v & v.wrapping_neg();
+        let r = v + c;
+        let n = (r ^ v).wrapping_shr(2 + v.trailing_zeros());
+        self.value = r | n;
+        Some(v)
+    }
+}
diff --git a/ext/crates/fp/src/prime/mod.rs b/ext/crates/fp/src/prime/mod.rs
new file mode 100644
index 000000000..58623952f
--- /dev/null
+++ b/ext/crates/fp/src/prime/mod.rs
@@ -0,0 +1,570 @@
+use std::{
+    fmt::Debug,
+    ops::{
+        Add, AddAssign, Div, DivAssign, Mul, MulAssign, Rem, RemAssign, Shl, ShlAssign, Shr,
+        ShrAssign, Sub, SubAssign,
+    },
+};
+
+use serde::{de::Error, Deserialize, Deserializer, Serialize, Serializer};
+
+pub mod binomial;
+pub mod iter;
+
+pub use binomial::Binomial;
+
+pub mod primes {
+    pub use super::{ValidPrime, P2, P3, P5, P7};
+}
+
+pub const TWO: ValidPrime = ValidPrime::new(2);
+
+/// A trait that represents a prime number. There are currently two kinds of structs that implement
+/// this trait: static primes and `ValidPrime`, the dynamic prime.
+///
+/// The methods in this trait take a `self` receiver so that the dynamic prime `ValidPrime` can
+/// implement it. We could also have a `&self` receiver, but since `Prime` is a supertrait of
+/// `Copy`, the two are equivalent. Using `self` might also be useful in the future if we ever want
+/// to play with autoref specialization.
+///
+/// The fact that e.g. `P2::to_u32` is hardcoded to return 2 means that a condition like `p.to_u32()
+/// == 2` (or even better, just `p == 2`) will reduce to `true` at compile time, and allow the
+/// compiler to eliminate an entire branch, while also leaving that check in for when the prime is
+/// chosen at runtime.
+pub trait Prime:
+    Debug
+    + Copy
+    + PartialEq
+    + Eq
+    + PartialEq<u32>
+    + PartialOrd<u32>
+    + Add<u32, Output = u32>
+    + Sub<u32, Output = u32>
+    + Mul<u32, Output = u32>
+    + Div<u32, Output = u32>
+    + Rem<u32, Output = u32>
+    + Shl<u32, Output = u32>
+    + Shr<u32, Output = u32>
+    + Serialize
+    + for<'de> Deserialize<'de>
+{
+    fn as_i32(self) -> i32;
+    fn to_dyn(self) -> ValidPrime;
+
+    fn as_u32(self) -> u32 {
+        self.as_i32() as u32
+    }
+
+    fn as_usize(self) -> usize {
+        self.as_u32() as usize
+    }
+
+    /// Computes the product mod p. This takes care of overflow.
+    fn product(self, n1: u32, n2: u32) -> u32 {
+        ((n1 as u64) * (n2 as u64) % (self.as_u32() as u64)) as u32
+    }
+
+    fn pow(self, exp: u32) -> u32 {
+        self.as_u32().pow(exp)
+    }
+
+    fn pow_mod(self, mut b: u32, mut e: u32) -> u32 {
+        assert!(self.as_u32() > 0);
+        let mut result: u32 = 1;
+        while e > 0 {
+            if (e & 1) == 1 {
+                result = self.product(result, b);
+            }
+            b = self.product(b, b);
+            e >>= 1;
+        }
+        result
+    }
+}
+
+#[derive(Debug, Clone, PartialEq, Eq)]
+pub enum PrimeError {
+    NotAnInteger(std::num::ParseIntError),
+    InvalidPrime(u32),
+}
+
+impl std::fmt::Display for PrimeError {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        match self {
+            PrimeError::NotAnInteger(s) => write!(f, "Not an integer: {}", s),
+            PrimeError::InvalidPrime(p) => write!(f, "{} is not a valid prime", p),
+        }
+    }
+}
+
+macro_rules! def_prime_static {
+    ($pn:ident, $p:literal) => {
+        #[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
+        pub struct $pn;
+
+        impl Prime for $pn {
+            #[inline]
+            fn as_i32(self) -> i32 {
+                $p
+            }
+            #[inline]
+            fn to_dyn(self) -> ValidPrime {
+                ValidPrime::new($p)
+            }
+        }
+
+        impl Serialize for $pn {
+            fn serialize<S>(&self, serializer: S) -> Result<S::Ok, S::Error>
+            where
+                S: Serializer,
+            {
+                self.as_u32().serialize(serializer)
+            }
+        }
+
+        impl<'de> Deserialize<'de> for $pn {
+            fn deserialize<D>(deserializer: D) -> Result<Self, D::Error>
+            where
+                D: Deserializer<'de>,
+            {
+                let p: u32 = u32::deserialize(deserializer)?;
+                $pn::try_from(p).map_err(D::Error::custom)
+            }
+        }
+    };
+}
+
+def_prime_static!(P2, 2);
+def_prime_static!(P3, 3);
+def_prime_static!(P5, 5);
+def_prime_static!(P7, 7);
+
+macro_rules! impl_op_pn_u32 {
+    ($pn:ty, $trt:ident, $mth:ident, $trt_assign:ident, $mth_assign:ident, $operator:tt) => {
+        impl $trt<$pn> for u32 {
+            type Output = u32;
+
+            fn $mth(self, other: $pn) -> Self::Output {
+                self $operator other.as_u32()
+            }
+        }
+
+        impl $trt<u32> for $pn {
+            type Output = u32;
+
+            fn $mth(self, other: u32) -> Self::Output {
+                self.as_u32() $operator other
+            }
+        }
+
+        impl $trt<$pn> for $pn {
+            type Output = u32;
+
+            fn $mth(self, other: $pn) -> Self::Output {
+                self.as_u32() $operator other.as_u32()
+            }
+        }
+
+        impl $trt_assign<$pn> for u32 {
+            fn $mth_assign(&mut self, other: $pn) {
+                *self = *self $operator other;
+            }
+        }
+    };
+}
+
+macro_rules! impl_prime_ops {
+    ($pn:ty) => {
+        impl_op_pn_u32!($pn, Add, add, AddAssign, add_assign, +);
+        impl_op_pn_u32!($pn, Sub, sub, SubAssign, sub_assign, -);
+        impl_op_pn_u32!($pn, Mul, mul, MulAssign, mul_assign, *);
+        impl_op_pn_u32!($pn, Div, div, DivAssign, div_assign, /);
+        impl_op_pn_u32!($pn, Rem, rem, RemAssign, rem_assign, %);
+        impl_op_pn_u32!($pn, Shl, shl, ShlAssign, shl_assign, <<);
+        impl_op_pn_u32!($pn, Shr, shr, ShrAssign, shr_assign, >>);
+
+        impl PartialEq<u32> for $pn {
+            fn eq(&self, other: &u32) -> bool {
+                self.as_u32() == *other
+            }
+        }
+
+        impl PartialOrd<u32> for $pn {
+            fn partial_cmp(&self, other: &u32) -> Option<std::cmp::Ordering> {
+                self.as_u32().partial_cmp(other)
+            }
+        }
+
+        impl From<$pn> for u32 {
+            fn from(value: $pn) -> u32 {
+                value.as_u32()
+            }
+        }
+
+        impl std::fmt::Display for $pn {
+            fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result {
+                <u32 as std::fmt::Display>::fmt(&self.as_u32(), f)
+            }
+        }
+    };
+}
+
+impl_prime_ops!(P2);
+impl_prime_ops!(P3);
+impl_prime_ops!(P5);
+impl_prime_ops!(P7);
+
+macro_rules! impl_try_from {
+    // We need the type both as a type and as an expression.
+    ($pn:ty, $pn_ex:expr) => {
+        impl std::convert::TryFrom<u32> for $pn {
+            type Error = PrimeError;
+
+            fn try_from(p: u32) -> Result<Self, PrimeError> {
+                if $pn_ex == p {
+                    Ok($pn_ex)
+                } else {
+                    Err(PrimeError::InvalidPrime(p))
+                }
+            }
+        }
+    };
+}
+
+impl_try_from!(P2, P2);
+impl_try_from!(P3, P3);
+impl_try_from!(P5, P5);
+impl_try_from!(P7, P7);
+
+#[cfg(feature = "odd-primes")]
+mod validprime {
+    use std::str::FromStr;
+
+    use super::*;
+
+    #[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
+    pub struct ValidPrime {
+        p: u32,
+    }
+
+    pub const fn is_prime(p: u32) -> bool {
+        // (2..p).all(|k| p % k != 0), but make it const
+        let mut k = 2;
+        while k < p {
+            if p % k == 0 {
+                return false;
+            }
+            k += 1;
+        }
+        true
+    }
+
+    impl ValidPrime {
+        pub const fn new(p: u32) -> ValidPrime {
+            // We need the size restriction for a few reasons.
+            //
+            // First, we need `bit_length(p)` to be smaller than 64. Otherwise, shifting a u64 by 64
+            // bits is considered an overflow. We could special case these shifts to be setting to
+            // 0, but that doesn't seem worth it.
+            //
+            // More importantly, the existence of `Prime::as_i32` means that we need `p` to fit in
+            // an i32. We want this method because there are a few places in the codebase that
+            // use it. It might be possible to go and change all of those to use `as_u32` instead,
+            // but it doesn't seem worth it for now.
+            assert!(p < (1 << 31), "Tried to construct a prime larger than 2^31");
+            assert!(is_prime(p), "Tried to construct a composite dynamic prime");
+            ValidPrime { p }
+        }
+
+        pub const fn new_unchecked(p: u32) -> ValidPrime {
+            ValidPrime { p }
+        }
+    }
+
+    impl Prime for ValidPrime {
+        fn as_i32(self) -> i32 {
+            self.p as i32
+        }
+
+        fn to_dyn(self) -> ValidPrime {
+            self
+        }
+    }
+
+    impl_prime_ops!(ValidPrime);
+
+    impl TryFrom<u32> for ValidPrime {
+        type Error = PrimeError;
+
+        fn try_from(p: u32) -> Result<Self, PrimeError> {
+            if is_prime(p) {
+                Ok(ValidPrime { p })
+            } else {
+                Err(PrimeError::InvalidPrime(p))
+            }
+        }
+    }
+
+    impl FromStr for ValidPrime {
+        type Err = PrimeError;
+
+        fn from_str(s: &str) -> Result<Self, Self::Err> {
+            let p: u32 = s.parse().map_err(PrimeError::NotAnInteger)?;
+            ValidPrime::try_from(p)
+        }
+    }
+
+    impl Serialize for ValidPrime {
+        fn serialize<S>(&self, serializer: S) -> Result<S::Ok, S::Error>
+        where
+            S: Serializer,
+        {
+            self.as_u32().serialize(serializer)
+        }
+    }
+
+    impl<'de> Deserialize<'de> for ValidPrime {
+        fn deserialize<D>(deserializer: D) -> Result<Self, D::Error>
+        where
+            D: Deserializer<'de>,
+        {
+            let p: u32 = u32::deserialize(deserializer)?;
+            ValidPrime::try_from(p).map_err(D::Error::custom)
+        }
+    }
+}
+
+#[cfg(not(feature = "odd-primes"))]
+mod validprime {
+    use std::str::FromStr;
+
+    use super::PrimeError;
+
+    pub type ValidPrime = super::P2;
+
+    pub const fn is_prime(p: u32) -> bool {
+        p == 2
+    }
+
+    impl ValidPrime {
+        pub const fn new(_p: u32) -> ValidPrime {
+            // Disregard the argument, assume the prime is 2. This has the advantage of us being
+            // able to use the same tests independently of whether odd-primes is enabled or not.
+            //
+            // This is sound but can cause some problems for the user that could be hard to
+            // diagnose. Maybe use debug_assert! and fix the tests?
+            super::P2
+        }
+
+        pub const fn new_unchecked(_p: u32) -> ValidPrime {
+            super::P2
+        }
+    }
+
+    impl FromStr for ValidPrime {
+        type Err = PrimeError;
+
+        fn from_str(s: &str) -> Result<Self, Self::Err> {
+            let p = s.parse().map_err(PrimeError::NotAnInteger)?;
+            if p == 2 {
+                Ok(super::P2)
+            } else {
+                Err(PrimeError::InvalidPrime(p))
+            }
+        }
+    }
+}
+
+pub use validprime::{is_prime, ValidPrime};
+
+/// Compute b^e mod p. This is a const version of `Prime::pow_mod`.
+pub const fn power_mod(p: u32, mut b: u32, mut e: u32) -> u32 {
+    // We can't use Prime::product because const traits are still unstable
+    assert!(p > 0);
+    let mut result: u32 = 1;
+    while e > 0 {
+        if (e & 1) == 1 {
+            result = ((result as u64) * (b as u64) % (p as u64)) as u32;
+        }
+        b = (((b as u64) * (b as u64)) % (p as u64)) as u32;
+        e >>= 1;
+    }
+    result
+}
+
+/// Compute the base 2 log of a number, rounded down to the nearest integer.
+///
+/// # Example
+/// ```
+/// # use fp::prime::log2;
+/// assert_eq!(0, log2(0b1));
+/// assert_eq!(1, log2(0b10));
+/// assert_eq!(1, log2(0b11));
+/// assert_eq!(3, log2(0b1011));
+/// ```
+pub const fn log2(n: usize) -> usize {
+    std::mem::size_of::<usize>() * 8 - 1 - n.leading_zeros() as usize
+}
+
+/// Discrete log base p of n.
+pub fn logp<P: Prime>(p: P, mut n: u32) -> u32 {
+    let mut result = 0u32;
+    while n > 0 {
+        n /= p.as_u32();
+        result += 1;
+    }
+    result
+}
+
+/// Factor $n$ as $p^k m$. Returns $(k, m)$.
+pub fn factor_pk<P: Prime>(p: P, mut n: u32) -> (u32, u32) {
+    if n == 0 {
+        return (0, 0);
+    }
+    let mut k = 0;
+    while n % p.as_u32() == 0 {
+        n /= p.as_u32();
+        k += 1;
+    }
+    (k, n)
+}
+
+// Uses a the lookup table we initialized.
+pub fn inverse<P: Prime>(p: P, k: u32) -> u32 {
+    use crate::constants::{INVERSE_TABLE, MAX_PRIME, PRIME_TO_INDEX_MAP};
+    assert!(k > 0 && p > k);
+
+    if p <= MAX_PRIME as u32 {
+        // LLVM doesn't understand the inequality is transitive
+        unsafe { *INVERSE_TABLE[PRIME_TO_INDEX_MAP[p.as_usize()]].get_unchecked(k as usize) }
+    } else {
+        power_mod(p.as_u32(), k, p.as_u32() - 2)
+    }
+}
+
+#[inline(always)]
+pub fn minus_one_to_the_n<P: Prime>(p: P, i: i32) -> u32 {
+    if i % 2 == 0 {
+        1
+    } else {
+        p - 1
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::{binomial::Binomial, inverse, iter::BinomialIterator, Prime, ValidPrime};
+    use crate::{
+        constants::PRIMES,
+        prime::{is_prime, PrimeError},
+    };
+
+    #[test]
+    fn validprime_test() {
+        for p in (0..(1 << 16)).filter(|&p| is_prime(p)) {
+            assert_eq!(ValidPrime::new(p), p);
+        }
+    }
+
+    #[test]
+    fn validprime_invalid() {
+        assert_eq!(
+            ValidPrime::try_from(4).unwrap_err(),
+            PrimeError::InvalidPrime(4)
+        );
+        assert_eq!(
+            "4".parse::<ValidPrime>().unwrap_err(),
+            PrimeError::InvalidPrime(4)
+        );
+        assert_eq!(
+            "4.0".parse::<ValidPrime>().unwrap_err(),
+            PrimeError::NotAnInteger("4.0".parse::<u32>().unwrap_err())
+        );
+    }
+
+    #[test]
+    fn inverse_test() {
+        for &p in PRIMES.iter() {
+            let p = ValidPrime::new(p);
+            for k in 1..p.as_u32() {
+                assert_eq!((inverse(p, k) * k) % p, 1);
+            }
+        }
+    }
+
+    #[test]
+    fn binomial_test() {
+        let entries = [[2, 2, 1, 0], [2, 3, 1, 1], [3, 1090, 730, 1], [7, 3, 2, 3]];
+
+        for entry in &entries {
+            assert_eq!(
+                entry[3],
+                u32::binomial(ValidPrime::new(entry[0]), entry[1], entry[2])
+            );
+        }
+    }
+
+    #[test]
+    fn binomial_vs_monomial() {
+        for &p in &[2, 3, 5, 7, 11] {
+            let p = ValidPrime::new(p);
+            for l in 0..20 {
+                for m in 0..20 {
+                    assert_eq!(u32::binomial(p, l + m, m), u32::multinomial(p, &mut [l, m]))
+                }
+            }
+        }
+    }
+
+    fn binomial_full(n: u32, j: u32) -> u32 {
+        let mut res = 1;
+        for k in j + 1..=n {
+            res *= k;
+        }
+        for k in 1..=(n - j) {
+            res /= k;
+        }
+        res
+    }
+
+    #[test]
+    fn binomial_cmp() {
+        for n in 0..12 {
+            for j in 0..=n {
+                let ans = binomial_full(n, j);
+                for &p in &[2, 3, 5, 7, 11] {
+                    assert_eq!(
+                        u32::binomial(ValidPrime::new(p), n, j),
+                        ans % p,
+                        "{n} choose {j} mod {p}"
+                    );
+                }
+                assert_eq!(u32::binomial4(n, j), ans % 4, "{n} choose {j} mod 4");
+                // binomial4_rec is only called on large n. It does not handle the n = 0, 1 cases
+                // correctly.
+                if n > 1 {
+                    assert_eq!(
+                        u32::binomial4_rec(n, j),
+                        ans % 4,
+                        "{n} choose {j} mod 4 rec"
+                    );
+                }
+            }
+        }
+    }
+
+    #[test]
+    fn binomial_iterator() {
+        let mut iter = BinomialIterator::new(4);
+        assert_eq!(iter.next(), Some(0b1111));
+        assert_eq!(iter.next(), Some(0b10111));
+        assert_eq!(iter.next(), Some(0b11011));
+        assert_eq!(iter.next(), Some(0b11101));
+        assert_eq!(iter.next(), Some(0b11110));
+        assert_eq!(iter.next(), Some(0b100111));
+        assert_eq!(iter.next(), Some(0b101011));
+        assert_eq!(iter.next(), Some(0b101101));
+        assert_eq!(iter.next(), Some(0b101110));
+        assert_eq!(iter.next(), Some(0b110011));
+    }
+}
diff --git a/ext/crates/fp/src/vector/impl_fpvectorp.rs b/ext/crates/fp/src/vector/impl_fpvectorp.rs
index 36dac5429..91d2a0e2e 100644
--- a/ext/crates/fp/src/vector/impl_fpvectorp.rs
+++ b/ext/crates/fp/src/vector/impl_fpvectorp.rs
@@ -2,7 +2,7 @@ use itertools::Itertools;
 
 use crate::{
     limb::{self, Limb},
-    prime::ValidPrime,
+    prime::{Prime, ValidPrime},
     simd,
 };
 
@@ -11,24 +11,25 @@ use super::{
     iter::{FpVectorIterator, FpVectorNonZeroIteratorP},
 };
 
-impl<const P: u32> FpVectorP<P> {
-    pub fn new_(len: usize) -> Self {
-        let number_of_limbs = limb::number::<P>(len);
+impl<P: Prime> FpVectorP<P> {
+    pub fn new(p: P, len: usize) -> Self {
+        let number_of_limbs = limb::number(p, len);
         Self {
+            p,
             len,
             limbs: vec![0; number_of_limbs],
         }
     }
 
-    pub fn from_raw_parts(len: usize, limbs: Vec<Limb>) -> Self {
-        debug_assert_eq!(limbs.len(), limb::number::<P>(len));
-        Self { len, limbs }
+    pub fn from_raw_parts(p: P, len: usize, limbs: Vec<Limb>) -> Self {
+        debug_assert_eq!(limbs.len(), limb::number(p, len));
+        Self { p, len, limbs }
     }
 
-    pub fn new_with_capacity_(len: usize, capacity: usize) -> Self {
-        let mut limbs = Vec::with_capacity(limb::number::<P>(capacity));
-        limbs.resize(limb::number::<P>(len), 0);
-        Self { len, limbs }
+    pub fn new_with_capacity(p: P, len: usize, capacity: usize) -> Self {
+        let mut limbs = Vec::with_capacity(limb::number(p, capacity));
+        limbs.resize(limb::number(p, len), 0);
+        Self { p, len, limbs }
     }
 
     pub const fn len(&self) -> usize {
@@ -39,14 +40,15 @@ impl<const P: u32> FpVectorP<P> {
         self.len == 0
     }
 
-    pub const fn prime(&self) -> ValidPrime {
-        ValidPrime::new(P)
+    pub fn prime(&self) -> ValidPrime {
+        self.p.to_dyn()
     }
 
     #[must_use]
     pub fn slice(&self, start: usize, end: usize) -> SliceP<'_, P> {
         assert!(start <= end && end <= self.len);
         SliceP {
+            p: self.p,
             limbs: &self.limbs,
             start,
             end,
@@ -57,6 +59,7 @@ impl<const P: u32> FpVectorP<P> {
     pub fn slice_mut(&mut self, start: usize, end: usize) -> SliceMutP<'_, P> {
         assert!(start <= end && end <= self.len);
         SliceMutP {
+            p: self.p,
             limbs: &mut self.limbs,
             start,
             end,
@@ -102,7 +105,7 @@ impl<const P: u32> FpVectorP<P> {
     }
 
     pub fn scale(&mut self, c: u32) {
-        match P {
+        match self.p.as_u32() {
             2 => {
                 if c == 0 {
                     self.set_to_zero()
@@ -110,12 +113,16 @@ impl<const P: u32> FpVectorP<P> {
             }
             3 | 5 => {
                 for limb in &mut self.limbs {
-                    *limb = limb::reduce::<P>(*limb * c as Limb);
+                    *limb = limb::reduce(self.p, *limb * c as Limb);
                 }
             }
             _ => {
                 for limb in &mut self.limbs {
-                    *limb = limb::pack::<_, P>(limb::unpack::<P>(*limb).map(|x| (x * c) % P));
+                    // We can cast x to u32 because we assume the limbs are reduced, so x < p < 2^31.
+                    *limb = limb::pack(
+                        self.p,
+                        limb::unpack(self.p, *limb).map(|x| self.p.product(x as u32, c)),
+                    );
                 }
             }
         }
@@ -125,17 +132,17 @@ impl<const P: u32> FpVectorP<P> {
     /// empty.
     pub fn add_offset(&mut self, other: &FpVectorP<P>, c: u32, offset: usize) {
         assert_eq!(self.len(), other.len());
-        let min_limb = offset / limb::entries_per_limb_const::<P>();
-        if P == 2 {
+        let min_limb = offset / limb::entries_per_limb(self.p);
+        if self.p == 2 {
             if c != 0 {
                 simd::add_simd(&mut self.limbs, &other.limbs, min_limb);
             }
         } else {
             for (left, right) in self.limbs.iter_mut().zip_eq(&other.limbs).skip(min_limb) {
-                *left = limb::add::<P>(*left, *right, c);
+                *left = limb::add(self.p, *left, *right, c);
             }
             for limb in &mut self.limbs[min_limb..] {
-                *limb = limb::reduce::<P>(*limb);
+                *limb = limb::reduce(self.p, *limb);
             }
         }
     }
@@ -144,8 +151,8 @@ impl<const P: u32> FpVectorP<P> {
     /// empty.
     pub fn add_offset_nosimd(&mut self, other: &FpVectorP<P>, c: u32, offset: usize) {
         assert_eq!(self.len(), other.len());
-        let min_limb = offset / limb::entries_per_limb_const::<P>();
-        if P == 2 {
+        let min_limb = offset / limb::entries_per_limb(self.p);
+        if self.p == 2 {
             if c != 0 {
                 for i in 0..self.limbs.len() {
                     self.limbs[i] ^= other.limbs[i];
@@ -153,10 +160,10 @@ impl<const P: u32> FpVectorP<P> {
             }
         } else {
             for (left, right) in self.limbs.iter_mut().zip_eq(&other.limbs).skip(min_limb) {
-                *left = limb::add::<P>(*left, *right, c);
+                *left = limb::add(self.p, *left, *right, c);
             }
             for limb in &mut self.limbs[min_limb..] {
-                *limb = limb::reduce::<P>(*limb);
+                *limb = limb::reduce(self.p, *limb);
             }
         }
     }
@@ -202,14 +209,14 @@ impl<const P: u32> FpVectorP<P> {
             return;
         }
         self.len = len;
-        self.limbs.resize(limb::number::<P>(len), 0);
+        self.limbs.resize(limb::number(self.p, len), 0);
     }
 
     /// This clears the vector and sets the length to `len`. This is useful for reusing
     /// allocations of temporary vectors.
     pub fn set_scratch_vector_size(&mut self, len: usize) {
         self.limbs.clear();
-        self.limbs.resize(limb::number::<P>(len), 0);
+        self.limbs.resize(limb::number(self.p, len), 0);
         self.len = len;
     }
 
@@ -221,8 +228,8 @@ impl<const P: u32> FpVectorP<P> {
         self.limbs.clear();
         self.limbs.extend(
             slice
-                .chunks(limb::entries_per_limb_const::<P>())
-                .map(|x| limb::pack::<_, P>(x.iter().copied())),
+                .chunks(limb::entries_per_limb(self.p))
+                .map(|x| limb::pack(self.p, x.iter().copied())),
         );
     }
 
@@ -230,7 +237,7 @@ impl<const P: u32> FpVectorP<P> {
     /// the number of entries per limb
     pub(crate) fn trim_start(&mut self, n: usize) {
         assert!(n <= self.len);
-        let entries_per = limb::entries_per_limb_const::<P>();
+        let entries_per = limb::entries_per_limb(self.p);
         assert_eq!(n % entries_per, 0);
         let num_limbs = n / entries_per;
         self.limbs.drain(0..num_limbs);
@@ -238,7 +245,7 @@ impl<const P: u32> FpVectorP<P> {
     }
 
     pub fn sign_rule(&self, other: &Self) -> bool {
-        assert_eq!(P, 2);
+        assert_eq!(self.p, 2);
         let mut result = 0;
         for target_limb_idx in 0..self.limbs.len() {
             let target_limb = other.limbs[target_limb_idx];
@@ -256,8 +263,8 @@ impl<const P: u32> FpVectorP<P> {
 
     pub fn add_truncate(&mut self, other: &Self, c: u32) -> Option<()> {
         for (left, right) in self.limbs.iter_mut().zip_eq(&other.limbs) {
-            *left = limb::add::<P>(*left, *right, c);
-            *left = limb::truncate::<P>(*left)?;
+            *left = limb::add(self.p, *left, *right, c);
+            *left = limb::truncate(self.p, *left)?;
         }
         Some(())
     }
@@ -266,7 +273,7 @@ impl<const P: u32> FpVectorP<P> {
     where
         for<'a> &'a mut T: TryInto<&'a mut Self>,
     {
-        if P == 2 {
+        if self.p == 2 {
             if c == 0 {
                 return false;
             }
@@ -306,9 +313,9 @@ impl<const P: u32> FpVectorP<P> {
 
     /// Find the index and value of the first non-zero entry of the vector. `None` if the vector is zero.
     pub fn first_nonzero(&self) -> Option<(usize, u32)> {
-        let entries_per_limb = limb::entries_per_limb_const::<P>();
-        let bit_length = limb::bit_length_const::<P>();
-        let bitmask = limb::bitmask::<P>();
+        let entries_per_limb = limb::entries_per_limb(self.p);
+        let bit_length = limb::bit_length(self.p);
+        let bitmask = limb::bitmask(self.p);
         for (i, &limb) in self.limbs.iter().enumerate() {
             if limb == 0 {
                 continue;
@@ -323,7 +330,7 @@ impl<const P: u32> FpVectorP<P> {
     }
 
     pub fn density(&self) -> f32 {
-        let num_nonzero = if P == 2 {
+        let num_nonzero = if self.p == 2 {
             self.limbs
                 .iter()
                 .copied()
@@ -336,21 +343,22 @@ impl<const P: u32> FpVectorP<P> {
     }
 }
 
-impl<T: AsRef<[u32]>, const P: u32> From<&T> for FpVectorP<P> {
-    fn from(slice: &T) -> Self {
-        let mut v = Self::new_(slice.as_ref().len());
+impl<T: AsRef<[u32]>, P: Prime> From<(P, &T)> for FpVectorP<P> {
+    fn from(data: (P, &T)) -> Self {
+        let (p, slice) = data;
+        let mut v = Self::new(p, slice.as_ref().len());
         v.limbs.clear();
         v.limbs.extend(
             slice
                 .as_ref()
-                .chunks(limb::entries_per_limb_const::<P>())
-                .map(|x| limb::pack::<_, P>(x.iter().copied())),
+                .chunks(limb::entries_per_limb(p))
+                .map(|x| limb::pack(p, x.iter().copied())),
         );
         v
     }
 }
 
-impl<const P: u32> From<&FpVectorP<P>> for Vec<u32> {
+impl<P: Prime> From<&FpVectorP<P>> for Vec<u32> {
     fn from(vec: &FpVectorP<P>) -> Vec<u32> {
         vec.iter().collect()
     }
diff --git a/ext/crates/fp/src/vector/impl_slicemutp.rs b/ext/crates/fp/src/vector/impl_slicemutp.rs
index cef7b040a..740d42d72 100644
--- a/ext/crates/fp/src/vector/impl_slicemutp.rs
+++ b/ext/crates/fp/src/vector/impl_slicemutp.rs
@@ -5,33 +5,33 @@ use itertools::Itertools;
 use crate::{
     constants,
     limb::{self, Limb},
-    prime::ValidPrime,
+    prime::{Prime, ValidPrime, P2},
 };
 
 use super::inner::{FpVectorP, SliceMutP, SliceP};
 
-impl<'a, const P: u32> SliceMutP<'a, P> {
+impl<'a, P: Prime> SliceMutP<'a, P> {
     pub fn prime(&self) -> ValidPrime {
-        ValidPrime::new(P)
+        self.p.to_dyn()
     }
 
     pub fn add_basis_element(&mut self, index: usize, value: u32) {
-        if P == 2 {
+        if self.p == 2 {
             // Checking for value % 2 == 0 appears to be less performant
-            let pair = limb::limb_bit_index_pair::<2>(index + self.start);
+            let pair = limb::limb_bit_index_pair(P2, index + self.start);
             self.limbs[pair.limb] ^= (value as Limb % 2) << pair.bit_index;
         } else {
             let mut x = self.as_slice().entry(index);
             x += value;
-            x %= P;
+            x %= self.p.as_u32();
             self.set_entry(index, x);
         }
     }
 
     pub fn set_entry(&mut self, index: usize, value: u32) {
         debug_assert!(index < self.as_slice().len());
-        let bit_mask = limb::bitmask::<P>();
-        let limb_index = limb::limb_bit_index_pair::<P>(index + self.start);
+        let bit_mask = limb::bitmask(self.p);
+        let limb_index = limb::limb_bit_index_pair(self.p, index + self.start);
         let mut result = self.limbs[limb_index.limb];
         result &= !(bit_mask << limb_index.bit_index);
         result |= (value as Limb) << limb_index.bit_index;
@@ -39,17 +39,17 @@ impl<'a, const P: u32> SliceMutP<'a, P> {
     }
 
     fn reduce_limbs(&mut self) {
-        if P != 2 {
+        if self.p != 2 {
             let limb_range = self.as_slice().limb_range();
 
             for limb in &mut self.limbs[limb_range] {
-                *limb = limb::reduce::<P>(*limb);
+                *limb = limb::reduce(self.p, *limb);
             }
         }
     }
 
     pub fn scale(&mut self, c: u32) {
-        if P == 2 {
+        if self.p == 2 {
             if c == 0 {
                 self.set_to_zero();
             }
@@ -97,12 +97,12 @@ impl<'a, const P: u32> SliceMutP<'a, P> {
     }
 
     pub fn add(&mut self, other: SliceP<'_, P>, c: u32) {
-        debug_assert!(c < P);
+        debug_assert!(self.p > c);
         if self.as_slice().is_empty() {
             return;
         }
 
-        if P == 2 {
+        if self.p == 2 {
             if c != 0 {
                 match self.as_slice().offset().cmp(&other.offset()) {
                     Ordering::Equal => self.add_shift_none(other, 1),
@@ -125,7 +125,7 @@ impl<'a, const P: u32> SliceMutP<'a, P> {
         let right_dim = right.len();
 
         for (i, v) in left.iter_nonzero() {
-            let entry = (v * coeff) % *self.prime();
+            let entry = (v * coeff) % self.prime();
             self.slice_mut(offset + i * right_dim, offset + (i + 1) * right_dim)
                 .add(right, entry);
         }
@@ -169,12 +169,13 @@ impl<'a, const P: u32> SliceMutP<'a, P> {
 
         let (min_mask, max_mask) = other.limb_masks();
 
-        self.limbs[target_range.start] = limb::add::<P>(
+        self.limbs[target_range.start] = limb::add(
+            self.p,
             self.limbs[target_range.start],
             other.limbs[source_range.start] & min_mask,
             c,
         );
-        self.limbs[target_range.start] = limb::reduce::<P>(self.limbs[target_range.start]);
+        self.limbs[target_range.start] = limb::reduce(self.p, self.limbs[target_range.start]);
 
         let target_inner_range = self.as_slice().limb_range_inner();
         let source_inner_range = other.limb_range_inner();
@@ -183,18 +184,20 @@ impl<'a, const P: u32> SliceMutP<'a, P> {
                 .iter_mut()
                 .zip_eq(&other.limbs[source_inner_range])
             {
-                *left = limb::add::<P>(*left, *right, c);
-                *left = limb::reduce::<P>(*left);
+                *left = limb::add(self.p, *left, *right, c);
+                *left = limb::reduce(self.p, *left);
             }
         }
         if source_range.len() > 1 {
             // The first and last limbs are distinct, so we process the last.
-            self.limbs[target_range.end - 1] = limb::add::<P>(
+            self.limbs[target_range.end - 1] = limb::add(
+                self.p,
                 self.limbs[target_range.end - 1],
                 other.limbs[source_range.end - 1] & max_mask,
                 c,
             );
-            self.limbs[target_range.end - 1] = limb::reduce::<P>(self.limbs[target_range.end - 1]);
+            self.limbs[target_range.end - 1] =
+                limb::reduce(self.p, self.limbs[target_range.end - 1]);
         }
     }
 
@@ -212,7 +215,7 @@ impl<'a, const P: u32> SliceMutP<'a, P> {
         }
 
         impl AddShiftLeftData {
-            fn new<const P: u32>(target: SliceP<'_, P>, source: SliceP<'_, P>) -> Self {
+            fn new<P: Prime>(p: P, target: SliceP<'_, P>, source: SliceP<'_, P>) -> Self {
                 debug_assert!(target.prime() == source.prime());
                 debug_assert!(target.offset() <= source.offset());
                 debug_assert!(
@@ -222,8 +225,8 @@ impl<'a, const P: u32> SliceMutP<'a, P> {
                     source.len()
                 );
                 let offset_shift = source.offset() - target.offset();
-                let bit_length = limb::bit_length_const::<P>();
-                let entries_per_limb = limb::entries_per_limb_const::<P>();
+                let bit_length = limb::bit_length(p);
+                let entries_per_limb = limb::entries_per_limb(p);
                 let usable_bits_per_limb = bit_length * entries_per_limb;
                 let tail_shift = usable_bits_per_limb - offset_shift;
                 let zero_bits = constants::BITS_PER_LIMB - usable_bits_per_limb;
@@ -248,73 +251,78 @@ impl<'a, const P: u32> SliceMutP<'a, P> {
                 }
             }
 
-            fn mask_first_limb<const P: u32>(&self, other: SliceP<'_, P>, i: usize) -> Limb {
+            fn mask_first_limb<P: Prime>(&self, other: SliceP<'_, P>, i: usize) -> Limb {
                 (other.limbs[i] & self.min_mask) >> self.offset_shift
             }
 
-            fn mask_middle_limb_a<const P: u32>(&self, other: SliceP<'_, P>, i: usize) -> Limb {
+            fn mask_middle_limb_a<P: Prime>(&self, other: SliceP<'_, P>, i: usize) -> Limb {
                 other.limbs[i] >> self.offset_shift
             }
 
-            fn mask_middle_limb_b<const P: u32>(&self, other: SliceP<'_, P>, i: usize) -> Limb {
+            fn mask_middle_limb_b<P: Prime>(&self, other: SliceP<'_, P>, i: usize) -> Limb {
                 (other.limbs[i] << (self.tail_shift + self.zero_bits)) >> self.zero_bits
             }
 
-            fn mask_last_limb_a<const P: u32>(&self, other: SliceP<'_, P>, i: usize) -> Limb {
+            fn mask_last_limb_a<P: Prime>(&self, other: SliceP<'_, P>, i: usize) -> Limb {
                 let source_limb_masked = other.limbs[i] & self.max_mask;
                 source_limb_masked << self.tail_shift
             }
 
-            fn mask_last_limb_b<const P: u32>(&self, other: SliceP<'_, P>, i: usize) -> Limb {
+            fn mask_last_limb_b<P: Prime>(&self, other: SliceP<'_, P>, i: usize) -> Limb {
                 let source_limb_masked = other.limbs[i] & self.max_mask;
                 source_limb_masked >> self.offset_shift
             }
         }
 
-        let dat = AddShiftLeftData::new(self.as_slice(), other);
+        let dat = AddShiftLeftData::new(self.p, self.as_slice(), other);
         let mut i = 0;
         {
-            self.limbs[i + dat.min_target_limb] = limb::add::<P>(
+            self.limbs[i + dat.min_target_limb] = limb::add(
+                self.p,
                 self.limbs[i + dat.min_target_limb],
                 dat.mask_first_limb(other, i + dat.min_source_limb),
                 c,
             );
         }
         for i in 1..dat.number_of_source_limbs - 1 {
-            self.limbs[i + dat.min_target_limb] = limb::add::<P>(
+            self.limbs[i + dat.min_target_limb] = limb::add(
+                self.p,
                 self.limbs[i + dat.min_target_limb],
                 dat.mask_middle_limb_a(other, i + dat.min_source_limb),
                 c,
             );
-            self.limbs[i + dat.min_target_limb - 1] = limb::add::<P>(
+            self.limbs[i + dat.min_target_limb - 1] = limb::add(
+                self.p,
                 self.limbs[i + dat.min_target_limb - 1],
                 dat.mask_middle_limb_b(other, i + dat.min_source_limb),
                 c,
             );
             self.limbs[i + dat.min_target_limb - 1] =
-                limb::reduce::<P>(self.limbs[i + dat.min_target_limb - 1]);
+                limb::reduce(self.p, self.limbs[i + dat.min_target_limb - 1]);
         }
         i = dat.number_of_source_limbs - 1;
         if i > 0 {
-            self.limbs[i + dat.min_target_limb - 1] = limb::add::<P>(
+            self.limbs[i + dat.min_target_limb - 1] = limb::add(
+                self.p,
                 self.limbs[i + dat.min_target_limb - 1],
                 dat.mask_last_limb_a(other, i + dat.min_source_limb),
                 c,
             );
             self.limbs[i + dat.min_target_limb - 1] =
-                limb::reduce::<P>(self.limbs[i + dat.min_target_limb - 1]);
+                limb::reduce(self.p, self.limbs[i + dat.min_target_limb - 1]);
             if dat.number_of_source_limbs == dat.number_of_target_limbs {
-                self.limbs[i + dat.min_target_limb] = limb::add::<P>(
+                self.limbs[i + dat.min_target_limb] = limb::add(
+                    self.p,
                     self.limbs[i + dat.min_target_limb],
                     dat.mask_last_limb_b(other, i + dat.min_source_limb),
                     c,
                 );
                 self.limbs[i + dat.min_target_limb] =
-                    limb::reduce::<P>(self.limbs[i + dat.min_target_limb]);
+                    limb::reduce(self.p, self.limbs[i + dat.min_target_limb]);
             }
         } else {
             self.limbs[i + dat.min_target_limb] =
-                limb::reduce::<P>(self.limbs[i + dat.min_target_limb]);
+                limb::reduce(self.p, self.limbs[i + dat.min_target_limb]);
         }
     }
 
@@ -332,7 +340,7 @@ impl<'a, const P: u32> SliceMutP<'a, P> {
         }
 
         impl AddShiftRightData {
-            fn new<const P: u32>(target: SliceP<'_, P>, source: SliceP<'_, P>) -> Self {
+            fn new<P: Prime>(p: P, target: SliceP<'_, P>, source: SliceP<'_, P>) -> Self {
                 debug_assert!(target.prime() == source.prime());
                 debug_assert!(target.offset() >= source.offset());
                 debug_assert!(
@@ -342,8 +350,8 @@ impl<'a, const P: u32> SliceMutP<'a, P> {
                     source.len()
                 );
                 let offset_shift = target.offset() - source.offset();
-                let bit_length = limb::bit_length_const::<P>();
-                let entries_per_limb = limb::entries_per_limb_const::<P>();
+                let bit_length = limb::bit_length(p);
+                let entries_per_limb = limb::entries_per_limb(p);
                 let usable_bits_per_limb = bit_length * entries_per_limb;
                 let tail_shift = usable_bits_per_limb - offset_shift;
                 let zero_bits = constants::BITS_PER_LIMB - usable_bits_per_limb;
@@ -367,47 +375,49 @@ impl<'a, const P: u32> SliceMutP<'a, P> {
                 }
             }
 
-            fn mask_first_limb_a<const P: u32>(&self, other: SliceP<'_, P>, i: usize) -> Limb {
+            fn mask_first_limb_a<P: Prime>(&self, other: SliceP<'_, P>, i: usize) -> Limb {
                 let source_limb_masked = other.limbs[i] & self.min_mask;
                 (source_limb_masked << (self.offset_shift + self.zero_bits)) >> self.zero_bits
             }
 
-            fn mask_first_limb_b<const P: u32>(&self, other: SliceP<'_, P>, i: usize) -> Limb {
+            fn mask_first_limb_b<P: Prime>(&self, other: SliceP<'_, P>, i: usize) -> Limb {
                 let source_limb_masked = other.limbs[i] & self.min_mask;
                 source_limb_masked >> self.tail_shift
             }
 
-            fn mask_middle_limb_a<const P: u32>(&self, other: SliceP<'_, P>, i: usize) -> Limb {
+            fn mask_middle_limb_a<P: Prime>(&self, other: SliceP<'_, P>, i: usize) -> Limb {
                 (other.limbs[i] << (self.offset_shift + self.zero_bits)) >> self.zero_bits
             }
 
-            fn mask_middle_limb_b<const P: u32>(&self, other: SliceP<'_, P>, i: usize) -> Limb {
+            fn mask_middle_limb_b<P: Prime>(&self, other: SliceP<'_, P>, i: usize) -> Limb {
                 other.limbs[i] >> self.tail_shift
             }
 
-            fn mask_last_limb_a<const P: u32>(&self, other: SliceP<'_, P>, i: usize) -> Limb {
+            fn mask_last_limb_a<P: Prime>(&self, other: SliceP<'_, P>, i: usize) -> Limb {
                 let source_limb_masked = other.limbs[i] & self.max_mask;
                 source_limb_masked << self.offset_shift
             }
 
-            fn mask_last_limb_b<const P: u32>(&self, other: SliceP<'_, P>, i: usize) -> Limb {
+            fn mask_last_limb_b<P: Prime>(&self, other: SliceP<'_, P>, i: usize) -> Limb {
                 let source_limb_masked = other.limbs[i] & self.max_mask;
                 source_limb_masked >> self.tail_shift
             }
         }
 
-        let dat = AddShiftRightData::new(self.as_slice(), other);
+        let dat = AddShiftRightData::new(self.p, self.as_slice(), other);
         let mut i = 0;
         {
-            self.limbs[i + dat.min_target_limb] = limb::add::<P>(
+            self.limbs[i + dat.min_target_limb] = limb::add(
+                self.p,
                 self.limbs[i + dat.min_target_limb],
                 dat.mask_first_limb_a(other, i + dat.min_source_limb),
                 c,
             );
             self.limbs[i + dat.min_target_limb] =
-                limb::reduce::<P>(self.limbs[i + dat.min_target_limb]);
+                limb::reduce(self.p, self.limbs[i + dat.min_target_limb]);
             if dat.number_of_target_limbs > 1 {
-                self.limbs[i + dat.min_target_limb + 1] = limb::add::<P>(
+                self.limbs[i + dat.min_target_limb + 1] = limb::add(
+                    self.p,
                     self.limbs[i + dat.min_target_limb + 1],
                     dat.mask_first_limb_b(other, i + dat.min_source_limb),
                     c,
@@ -415,14 +425,16 @@ impl<'a, const P: u32> SliceMutP<'a, P> {
             }
         }
         for i in 1..dat.number_of_source_limbs - 1 {
-            self.limbs[i + dat.min_target_limb] = limb::add::<P>(
+            self.limbs[i + dat.min_target_limb] = limb::add(
+                self.p,
                 self.limbs[i + dat.min_target_limb],
                 dat.mask_middle_limb_a(other, i + dat.min_source_limb),
                 c,
             );
             self.limbs[i + dat.min_target_limb] =
-                limb::reduce::<P>(self.limbs[i + dat.min_target_limb]);
-            self.limbs[i + dat.min_target_limb + 1] = limb::add::<P>(
+                limb::reduce(self.p, self.limbs[i + dat.min_target_limb]);
+            self.limbs[i + dat.min_target_limb + 1] = limb::add(
+                self.p,
                 self.limbs[i + dat.min_target_limb + 1],
                 dat.mask_middle_limb_b(other, i + dat.min_source_limb),
                 c,
@@ -430,15 +442,17 @@ impl<'a, const P: u32> SliceMutP<'a, P> {
         }
         i = dat.number_of_source_limbs - 1;
         if i > 0 {
-            self.limbs[i + dat.min_target_limb] = limb::add::<P>(
+            self.limbs[i + dat.min_target_limb] = limb::add(
+                self.p,
                 self.limbs[i + dat.min_target_limb],
                 dat.mask_last_limb_a(other, i + dat.min_source_limb),
                 c,
             );
             self.limbs[i + dat.min_target_limb] =
-                limb::reduce::<P>(self.limbs[i + dat.min_target_limb]);
+                limb::reduce(self.p, self.limbs[i + dat.min_target_limb]);
             if dat.number_of_target_limbs > dat.number_of_source_limbs {
-                self.limbs[i + dat.min_target_limb + 1] = limb::add::<P>(
+                self.limbs[i + dat.min_target_limb + 1] = limb::add(
+                    self.p,
                     self.limbs[i + dat.min_target_limb + 1],
                     dat.mask_last_limb_b(other, i + dat.min_source_limb),
                     c,
@@ -447,7 +461,7 @@ impl<'a, const P: u32> SliceMutP<'a, P> {
         }
         if dat.number_of_target_limbs > dat.number_of_source_limbs {
             self.limbs[i + dat.min_target_limb + 1] =
-                limb::reduce::<P>(self.limbs[i + dat.min_target_limb + 1]);
+                limb::reduce(self.p, self.limbs[i + dat.min_target_limb + 1]);
         }
     }
 
@@ -475,6 +489,7 @@ impl<'a, const P: u32> SliceMutP<'a, P> {
         assert!(start <= end && end <= self.as_slice().len());
 
         SliceMutP {
+            p: self.p,
             limbs: &mut *self.limbs,
             start: self.start + start,
             end: self.start + end,
@@ -485,6 +500,7 @@ impl<'a, const P: u32> SliceMutP<'a, P> {
     #[must_use]
     pub fn as_slice(&self) -> SliceP<'_, P> {
         SliceP {
+            p: self.p,
             limbs: &*self.limbs,
             start: self.start,
             end: self.end,
@@ -496,6 +512,7 @@ impl<'a, const P: u32> SliceMutP<'a, P> {
     #[must_use]
     pub fn copy(&mut self) -> SliceMutP<'_, P> {
         SliceMutP {
+            p: self.p,
             limbs: self.limbs,
             start: self.start,
             end: self.end,
@@ -503,7 +520,7 @@ impl<'a, const P: u32> SliceMutP<'a, P> {
     }
 }
 
-impl<'a, const P: u32> From<&'a mut FpVectorP<P>> for SliceMutP<'a, P> {
+impl<'a, P: Prime> From<&'a mut FpVectorP<P>> for SliceMutP<'a, P> {
     fn from(v: &'a mut FpVectorP<P>) -> Self {
         v.slice_mut(0, v.len)
     }
diff --git a/ext/crates/fp/src/vector/impl_slicep.rs b/ext/crates/fp/src/vector/impl_slicep.rs
index 3e2bc51df..fff8cd12d 100644
--- a/ext/crates/fp/src/vector/impl_slicep.rs
+++ b/ext/crates/fp/src/vector/impl_slicep.rs
@@ -1,7 +1,7 @@
 use crate::{
     constants,
     limb::{self, Limb},
-    prime::ValidPrime,
+    prime::{Prime, ValidPrime},
 };
 
 use super::{
@@ -11,9 +11,9 @@ use super::{
 
 // Public methods
 
-impl<'a, const P: u32> SliceP<'a, P> {
+impl<'a, P: Prime> SliceP<'a, P> {
     pub fn prime(&self) -> ValidPrime {
-        ValidPrime::new(P)
+        self.p.to_dyn()
     }
 
     pub fn len(&self) -> usize {
@@ -31,8 +31,8 @@ impl<'a, const P: u32> SliceP<'a, P> {
             index,
             self.len()
         );
-        let bit_mask = limb::bitmask::<P>();
-        let limb_index = limb::limb_bit_index_pair::<P>(index + self.start);
+        let bit_mask = limb::bitmask(self.p);
+        let limb_index = limb::limb_bit_index_pair(self.p, index + self.start);
         let mut result = self.limbs[limb_index.limb];
         result >>= limb_index.bit_index;
         result &= bit_mask;
@@ -72,6 +72,7 @@ impl<'a, const P: u32> SliceP<'a, P> {
         assert!(start <= end && end <= self.len());
 
         SliceP {
+            p: self.p,
             limbs: self.limbs,
             start: self.start + start,
             end: self.start + end,
@@ -81,8 +82,8 @@ impl<'a, const P: u32> SliceP<'a, P> {
     /// Converts a slice to an owned FpVectorP. This is vastly more efficient if the start of the vector is aligned.
     #[must_use]
     pub fn to_owned(self) -> FpVectorP<P> {
-        let mut new = FpVectorP::<P>::new_(self.len());
-        if self.start % limb::entries_per_limb_const::<P>() == 0 {
+        let mut new = FpVectorP::new(self.p, self.len());
+        if self.start % limb::entries_per_limb(self.p) == 0 {
             let limb_range = self.limb_range();
             new.limbs[0..limb_range.len()].copy_from_slice(&self.limbs[limb_range]);
             if !new.limbs.is_empty() {
@@ -97,17 +98,17 @@ impl<'a, const P: u32> SliceP<'a, P> {
 }
 
 // Limb methods
-impl<'a, const P: u32> SliceP<'a, P> {
+impl<'a, P: Prime> SliceP<'a, P> {
     #[inline]
     pub(super) fn offset(&self) -> usize {
-        let bit_length = limb::bit_length_const::<P>();
-        let entries_per_limb = limb::entries_per_limb_const::<P>();
+        let bit_length = limb::bit_length(self.p);
+        let entries_per_limb = limb::entries_per_limb(self.p);
         (self.start % entries_per_limb) * bit_length
     }
 
     #[inline]
     pub(super) fn limb_range(&self) -> std::ops::Range<usize> {
-        limb::range::<P>(self.start, self.end)
+        limb::range(self.p, self.start, self.end)
     }
 
     /// This function underflows if `self.end == 0`, which happens if and only if we are taking a
@@ -127,8 +128,8 @@ impl<'a, const P: u32> SliceP<'a, P> {
 
     #[inline(always)]
     pub(super) fn max_limb_mask(&self) -> Limb {
-        let num_entries = 1 + (self.end - 1) % limb::entries_per_limb_const::<P>();
-        let bit_max = num_entries * limb::bit_length_const::<P>();
+        let num_entries = 1 + (self.end - 1) % limb::entries_per_limb(self.p);
+        let bit_max = num_entries * limb::bit_length(self.p);
 
         (!0) >> (constants::BITS_PER_LIMB - bit_max)
     }
@@ -146,7 +147,7 @@ impl<'a, const P: u32> SliceP<'a, P> {
     }
 }
 
-impl<'a, const P: u32> From<&'a FpVectorP<P>> for SliceP<'a, P> {
+impl<'a, P: Prime> From<&'a FpVectorP<P>> for SliceP<'a, P> {
     fn from(v: &'a FpVectorP<P>) -> Self {
         v.slice(0, v.len)
     }
diff --git a/ext/crates/fp/src/vector/inner.rs b/ext/crates/fp/src/vector/inner.rs
index b47bd838a..3e4afd070 100644
--- a/ext/crates/fp/src/vector/inner.rs
+++ b/ext/crates/fp/src/vector/inner.rs
@@ -1,7 +1,7 @@
 // This generates better llvm optimization
 #![allow(clippy::int_plus_one)]
 
-use crate::limb::Limb;
+use crate::{limb::Limb, prime::Prime};
 
 /// An `FpVectorP` is a vector over $\mathbb{F}_p$ for a fixed prime, implemented using const
 /// generics. Due to limitations with const generics, we cannot constrain P to actually be a prime,
@@ -10,14 +10,16 @@ use crate::limb::Limb;
 /// Interally, it packs entries of the vectors into limbs. However, this is an abstraction that
 /// must not leave the `fp` library.
 #[derive(Debug, Hash, Eq, PartialEq, Clone)]
-pub struct FpVectorP<const P: u32> {
+pub struct FpVectorP<P: Prime> {
+    pub(super) p: P,
     pub(super) len: usize,
     pub(super) limbs: Vec<Limb>,
 }
 
 /// A SliceP is a slice of an FpVectorP. This immutably borrows the vector and implements Copy
 #[derive(Debug, Copy, Clone)]
-pub struct SliceP<'a, const P: u32> {
+pub struct SliceP<'a, P: Prime> {
+    pub(super) p: P,
     pub(super) limbs: &'a [Limb],
     pub(super) start: usize,
     pub(super) end: usize,
@@ -28,7 +30,8 @@ pub struct SliceP<'a, const P: u32> {
 /// that imitates the reborrowing, that mutably borrows `SliceMutP` and returns a `SliceMutP` with
 /// a shorter lifetime.
 #[derive(Debug)]
-pub struct SliceMutP<'a, const P: u32> {
+pub struct SliceMutP<'a, P: Prime> {
+    pub(super) p: P,
     pub(super) limbs: &'a mut [Limb],
     pub(super) start: usize,
     pub(super) end: usize,
diff --git a/ext/crates/fp/src/vector/iter.rs b/ext/crates/fp/src/vector/iter.rs
index 94bf88533..0e62abc05 100644
--- a/ext/crates/fp/src/vector/iter.rs
+++ b/ext/crates/fp/src/vector/iter.rs
@@ -1,4 +1,7 @@
-use crate::limb::{self, Limb};
+use crate::{
+    limb::{self, Limb},
+    prime::Prime,
+};
 
 use super::inner::SliceP;
 
@@ -14,7 +17,7 @@ pub struct FpVectorIterator<'a> {
 }
 
 impl<'a> FpVectorIterator<'a> {
-    pub(super) fn new<const P: u32>(vec: SliceP<'a, P>) -> Self {
+    pub(super) fn new<P: Prime>(vec: SliceP<'a, P>) -> Self {
         let counter = vec.len();
         let limbs = &vec.limbs;
 
@@ -30,17 +33,17 @@ impl<'a> FpVectorIterator<'a> {
                 counter,
             };
         }
-        let pair = limb::limb_bit_index_pair::<P>(vec.start);
+        let pair = limb::limb_bit_index_pair(vec.p, vec.start);
 
-        let bit_length = limb::bit_length_const::<P>();
+        let bit_length = limb::bit_length(vec.p);
         let cur_limb = limbs[pair.limb] >> pair.bit_index;
 
-        let entries_per_limb = limb::entries_per_limb_const::<P>();
+        let entries_per_limb = limb::entries_per_limb(vec.p);
         Self {
             limbs,
             bit_length,
             entries_per_limb_m_1: entries_per_limb - 1,
-            bit_mask: limb::bitmask::<P>(),
+            bit_mask: limb::bitmask(vec.p),
             limb_index: pair.limb,
             entries_left: entries_per_limb - (vec.start % entries_per_limb),
             cur_limb,
@@ -108,7 +111,8 @@ impl<'a> ExactSizeIterator for FpVectorIterator<'a> {
 
 /// Iterator over non-zero entries of an FpVector. This is monomorphized over P for significant
 /// performance gains.
-pub struct FpVectorNonZeroIteratorP<'a, const P: u32> {
+pub struct FpVectorNonZeroIteratorP<'a, P> {
+    p: P,
     limbs: &'a [Limb],
     limb_index: usize,
     cur_limb_entries_left: usize,
@@ -117,15 +121,16 @@ pub struct FpVectorNonZeroIteratorP<'a, const P: u32> {
     dim: usize,
 }
 
-impl<'a, const P: u32> FpVectorNonZeroIteratorP<'a, P> {
+impl<'a, P: Prime> FpVectorNonZeroIteratorP<'a, P> {
     pub(super) fn new(vec: SliceP<'a, P>) -> Self {
-        let entries_per_limb = limb::entries_per_limb_const::<P>();
+        let entries_per_limb = limb::entries_per_limb(vec.p);
 
         let dim = vec.len();
         let limbs = vec.limbs;
 
         if dim == 0 {
             return Self {
+                p: vec.p,
                 limbs,
                 limb_index: 0,
                 cur_limb_entries_left: 0,
@@ -135,10 +140,11 @@ impl<'a, const P: u32> FpVectorNonZeroIteratorP<'a, P> {
             };
         }
         let min_index = vec.start;
-        let pair = limb::limb_bit_index_pair::<P>(min_index);
+        let pair = limb::limb_bit_index_pair(vec.p, min_index);
         let cur_limb = limbs[pair.limb] >> pair.bit_index;
         let cur_limb_entries_left = entries_per_limb - (min_index % entries_per_limb);
         Self {
+            p: vec.p,
             limbs,
             limb_index: pair.limb,
             cur_limb_entries_left,
@@ -149,12 +155,12 @@ impl<'a, const P: u32> FpVectorNonZeroIteratorP<'a, P> {
     }
 }
 
-impl<'a, const P: u32> Iterator for FpVectorNonZeroIteratorP<'a, P> {
+impl<'a, P: Prime> Iterator for FpVectorNonZeroIteratorP<'a, P> {
     type Item = (usize, u32);
     fn next(&mut self) -> Option<Self::Item> {
-        let bit_length: usize = limb::bit_length_const::<P>();
-        let bitmask: Limb = limb::bitmask::<P>();
-        let entries_per_limb: usize = limb::entries_per_limb_const::<P>();
+        let bit_length: usize = limb::bit_length(self.p);
+        let bitmask: Limb = limb::bitmask(self.p);
+        let entries_per_limb: usize = limb::entries_per_limb(self.p);
         loop {
             let bits_left = (self.cur_limb_entries_left * bit_length) as u32;
             #[allow(clippy::unnecessary_cast)]
diff --git a/ext/crates/fp/src/vector/mod.rs b/ext/crates/fp/src/vector/mod.rs
index 8dbca234f..efdc9e601 100644
--- a/ext/crates/fp/src/vector/mod.rs
+++ b/ext/crates/fp/src/vector/mod.rs
@@ -17,13 +17,15 @@ mod iter;
 
 #[cfg(test)]
 mod test {
+    use std::sync::OnceLock;
+
     use itertools::Itertools;
     use proptest::prelude::*;
     use rstest::rstest;
 
     use super::{inner::FpVectorP, *};
-    use crate::limb;
-    use crate::prime::ValidPrime;
+    use crate::limb::{self, bit_length};
+    use crate::prime::{Prime, ValidPrime, P2};
 
     pub struct VectorDiffEntry {
         pub index: usize,
@@ -100,14 +102,30 @@ mod test {
         (0..dimension).map(|_| rng.gen_range(0..p)).collect()
     }
 
-    /// An arbitrary `ValidPrime`
+    /// An arbitrary `ValidPrime` in the range `2..(1 << 24)`, plus the largest prime that we support.
     fn arb_prime() -> impl Strategy<Value = ValidPrime> {
-        prop_oneof![
-            Just(ValidPrime::new(2)),
-            Just(ValidPrime::new(3)),
-            Just(ValidPrime::new(5)),
-            Just(ValidPrime::new(7)),
-        ]
+        static TEST_PRIMES: OnceLock<Vec<ValidPrime>> = OnceLock::new();
+        let test_primes = TEST_PRIMES.get_or_init(|| {
+            // Sieve of erathosthenes
+            const MAX: usize = 1 << 24;
+            let mut is_prime = Vec::new();
+            is_prime.resize_with(MAX, || true);
+            is_prime[0] = false;
+            is_prime[1] = false;
+            for i in 2..MAX {
+                if is_prime[i] {
+                    for j in ((2 * i)..MAX).step_by(i) {
+                        is_prime[j] = false;
+                    }
+                }
+            }
+            (0..MAX)
+                .filter(|&i| is_prime[i])
+                .map(|p| ValidPrime::new_unchecked(p as u32))
+                .chain(std::iter::once(ValidPrime::new_unchecked(2147483647)))
+                .collect()
+        });
+        (0..test_primes.len()).prop_map(|i| test_primes[i])
     }
 
     /// An arbitrary (prime, dimension) pair
@@ -135,7 +153,7 @@ mod test {
     /// take in a `Vec<u32>` instead of an `FpVector` directly because they will usually apply some
     /// operation on both the `FpVector` and the original `Vec<u32>` and then compare the results.
     fn arb_vec_u32(p: ValidPrime, dimension: usize) -> impl Strategy<Value = Vec<u32>> {
-        proptest::collection::vec(0..*p, dimension)
+        proptest::collection::vec(0..p.as_u32(), dimension)
     }
 
     /// A pair of a prime `p` and a vector containing values in the range `0..p`. In other
@@ -193,6 +211,11 @@ mod test {
             .. ProptestConfig::default()
         })]
 
+        #[test]
+        fn test_bit_length(p in arb_prime()) {
+            prop_assert!(bit_length(p) <= 63);
+        }
+
         #[test]
         fn test_serialize((p, v_arr) in arb_vec()) {
             use std::io::{Seek, Cursor};
@@ -215,19 +238,19 @@ mod test {
             v.add(&w, 1);
 
             for (v_element, w_element) in v_arr.iter_mut().zip(w_arr.iter()) {
-                *v_element = (*v_element + *w_element) % *p;
+                *v_element = (*v_element + *w_element) % p;
             }
             v.assert_list_eq(&v_arr);
         }
 
         #[test]
         fn test_scale((p, mut v_arr, c) in arb_prime_dim().prop_flat_map(|(p, dim)| {
-            (Just(p), arb_vec_u32(p, dim), 0u32..*p)
+            (Just(p), arb_vec_u32(p, dim), 0..p.as_u32())
         })) {
             let mut v = FpVector::from_slice(p, &v_arr);
             v.scale(c);
             for entry in &mut v_arr {
-                *entry = (*entry * c) % *p;
+                *entry = p.product(*entry, c);
             }
             v.assert_list_eq(&v_arr);
         }
@@ -235,14 +258,14 @@ mod test {
         #[test]
         fn test_scale_slice((p, mut v_arr, (slice_start, slice_end), c) in
             arb_prime_dim().prop_flat_map(|(p, dim)| {
-                (Just(p), arb_vec_u32(p, dim), arb_slice(dim), 0u32..*p)
+                (Just(p), arb_vec_u32(p, dim), arb_slice(dim), 0..p.as_u32())
             })
         ) {
             let mut v = FpVector::from_slice(p, &v_arr);
             v.slice_mut(slice_start, slice_end).scale(c);
 
             for entry in &mut v_arr[slice_start..slice_end] {
-                *entry = (*entry * c) % *p;
+                *entry = p.product(*entry, c);
             }
             v.assert_list_eq(&v_arr);
         }
@@ -331,7 +354,7 @@ mod test {
                 .add(w.slice(slice_start, slice_end), 1);
 
             for i in slice_start..slice_end {
-                v_arr[i] = (v_arr[i] + w_arr[i]) % *p;
+                v_arr[i] = (v_arr[i] + w_arr[i]) % p;
             }
             v.assert_list_eq(&v_arr);
         }
@@ -389,7 +412,7 @@ mod test {
                 .iter_mut()
                 .zip(w_arr[slice2_start..slice2_end].iter())
             {
-                *v_element = (*v_element + *w_element) % *p;
+                *v_element = (*v_element + *w_element) % p;
             }
             v.assert_list_eq(&v_arr);
         }
@@ -403,7 +426,7 @@ mod test {
 
             for (i, x) in v_small.iter_mut().enumerate() {
                 *x += v_big[mask[i]];
-                *x %= *p;
+                *x %= p;
             }
 
             v.assert_list_eq(&v_small);
@@ -417,7 +440,7 @@ mod test {
             v.as_slice_mut().add_unmasked(w.as_slice(), 1, &mask);
             for (i, &x) in v_small.iter().enumerate() {
                 v_big[mask[i]] += x;
-                v_big[mask[i]] %= *p;
+                v_big[mask[i]] %= p;
             }
             v.assert_list_eq(&v_big);
         }
@@ -527,7 +550,7 @@ mod test {
 
         let p = ValidPrime::new(p);
         const E_MAX: usize = 4;
-        let pto_the_e_max = (*p * *p * *p * *p) * *p;
+        let pto_the_e_max = (p * p * p * p) * p;
         let mut v = Vec::with_capacity(E_MAX + 1);
         let mut w = Vec::with_capacity(E_MAX + 1);
         for _ in 0..=E_MAX {
@@ -540,10 +563,10 @@ mod test {
             let mut ev = v_arr[i];
             let mut ew = w_arr[i];
             for e in 0..=E_MAX {
-                v[e].set_entry(i, ev % *p);
-                w[e].set_entry(i, ew % *p);
-                ev /= *p;
-                ew /= *p;
+                v[e].set_entry(i, ev % p);
+                w[e].set_entry(i, ew % p);
+                ev /= p;
+                ew /= p;
             }
         }
 
@@ -567,7 +590,7 @@ mod test {
         let mut vec_result = vec![0; dim];
         for (i, entry) in vec_result.iter_mut().enumerate() {
             for e in (0..=E_MAX).rev() {
-                *entry *= *p;
+                *entry *= p;
                 *entry += v[e].entry(i);
             }
         }
@@ -609,8 +632,8 @@ mod test {
     #[test]
     #[ignore]
     fn test_sign_rule() {
-        let mut in1 = FpVectorP::<2>::new_(128);
-        let mut in2 = FpVectorP::<2>::new_(128);
+        let mut in1 = FpVectorP::new(P2, 128);
+        let mut in2 = FpVectorP::new(P2, 128);
         let tests = [
             (
                 0x181e20846a820820,
diff --git a/ext/crates/fp/src/vector/vector_2.rs b/ext/crates/fp/src/vector/vector_2.rs
index 40f89ba03..c56effc06 100644
--- a/ext/crates/fp/src/vector/vector_2.rs
+++ b/ext/crates/fp/src/vector/vector_2.rs
@@ -8,40 +8,36 @@ use itertools::Itertools;
 
 use serde::{Deserialize, Deserializer, Serialize, Serializer};
 
-use crate::limb::{entries_per_limb_const, Limb};
-use crate::prime::ValidPrime;
+use crate::limb::{entries_per_limb, Limb};
+use crate::prime::{Prime, ValidPrime, P2};
 use crate::vector::inner::{FpVectorP, SliceMutP, SliceP};
 
 use super::iter::{FpVectorIterator, FpVectorNonZeroIteratorP};
 
-pub type FpVector = FpVectorP<2>;
-pub type Slice<'a> = SliceP<'a, 2>;
-pub type SliceMut<'a> = SliceMutP<'a, 2>;
-pub type FpVectorNonZeroIterator<'a> = FpVectorNonZeroIteratorP<'a, 2>;
+pub type FpVector = FpVectorP<P2>;
+pub type Slice<'a> = SliceP<'a, P2>;
+pub type SliceMut<'a> = SliceMutP<'a, P2>;
+pub type FpVectorNonZeroIterator<'a> = FpVectorNonZeroIteratorP<'a, P2>;
 
 // `FpVector` implementations
 
 impl FpVector {
-    pub fn new(_p: ValidPrime, len: usize) -> FpVector {
-        FpVector::new_(len)
-    }
-
-    pub fn new_with_capacity(_p: ValidPrime, len: usize, capacity: usize) -> FpVector {
-        FpVector::new_with_capacity_(len, capacity)
-    }
-
-    pub fn from_slice(_p: ValidPrime, slice: &[u32]) -> Self {
-        Self::from(&slice)
+    pub fn from_slice<P: Prime>(p: P, slice: &[u32]) -> Self {
+        if p == 2 {
+            Self::from((P2, &slice))
+        } else {
+            panic!("Only p = 2 is supported")
+        }
     }
 
     pub fn num_limbs(_p: ValidPrime, len: usize) -> usize {
-        let entries_per_limb = entries_per_limb_const::<2>();
+        let entries_per_limb = entries_per_limb(P2);
         (len + entries_per_limb - 1) / entries_per_limb
     }
 
     #[allow(dead_code)]
     pub(crate) fn padded_len(p: ValidPrime, len: usize) -> usize {
-        Self::num_limbs(p, len) * entries_per_limb_const::<2>()
+        Self::num_limbs(p, len) * entries_per_limb(P2)
     }
 
     pub fn update_from_bytes(&mut self, data: &mut impl Read) -> std::io::Result<()> {
diff --git a/ext/crates/fp/src/vector/vector_generic.rs b/ext/crates/fp/src/vector/vector_generic.rs
index 647988d54..94661d6c5 100644
--- a/ext/crates/fp/src/vector/vector_generic.rs
+++ b/ext/crates/fp/src/vector/vector_generic.rs
@@ -16,7 +16,7 @@ use itertools::Itertools;
 use serde::{Deserialize, Deserializer, Serialize, Serializer};
 
 use crate::limb::{entries_per_limb, Limb};
-use crate::prime::ValidPrime;
+use crate::prime::{primes::*, Prime, ValidPrime};
 use crate::vector::inner::{FpVectorP, SliceMutP, SliceP};
 
 use super::iter::{FpVectorIterator, FpVectorNonZeroIteratorP};
@@ -31,6 +31,7 @@ macro_rules! dispatch_vector_inner {
                 (Self::_3(x), $other::_3(y)) => x.$method(y, $($arg),*),
                 (Self::_5(x), $other::_5(y)) => x.$method(y, $($arg),*),
                 (Self::_7(x), $other::_7(y)) => x.$method(y, $($arg),*),
+                (Self::Big(x), $other::Big(y)) => x.$method(y, $($arg),*),
                 (l, r) => {
                     panic!("Applying {} to vectors over different primes ({} and {})", stringify!($method), l.prime(), r.prime());
                 }
@@ -45,6 +46,7 @@ macro_rules! dispatch_vector_inner {
                 (Self::_3(x), $other::_3(y)) => x.$method(y, $($arg),*),
                 (Self::_5(x), $other::_5(y)) => x.$method(y, $($arg),*),
                 (Self::_7(x), $other::_7(y)) => x.$method(y, $($arg),*),
+                (Self::Big(x), $other::Big(y)) => x.$method(y, $($arg),*),
                 (l, r) => {
                     panic!("Applying {} to vectors over different primes ({} and {})", stringify!($method), l.prime(), r.prime());
                 }
@@ -58,6 +60,7 @@ macro_rules! dispatch_vector_inner {
                 (Self::_3(x), $other::_3(y)) => x.$method(y, $($arg),*),
                 (Self::_5(x), $other::_5(y)) => x.$method(y, $($arg),*),
                 (Self::_7(x), $other::_7(y)) => x.$method(y, $($arg),*),
+                (Self::Big(x), $other::Big(y)) => x.$method(y, $($arg),*),
                 (l, r) => {
                     panic!("Applying {} to vectors over different primes ({} and {})", stringify!($method), l.prime(), r.prime());
                 }
@@ -72,6 +75,7 @@ macro_rules! dispatch_vector_inner {
                 Self::_3(x) => $ret::_3(x.$method($($arg),*)),
                 Self::_5(x) => $ret::_5(x.$method($($arg),*)),
                 Self::_7(x) => $ret::_7(x.$method($($arg),*)),
+                Self::Big(x) => $ret::Big(x.$method($($arg),*)),
             }
         }
     };
@@ -83,6 +87,7 @@ macro_rules! dispatch_vector_inner {
                 Self::_3(x) => $ret::_3(x.$method($($arg),*)),
                 Self::_5(x) => $ret::_5(x.$method($($arg),*)),
                 Self::_7(x) => $ret::_7(x.$method($($arg),*)),
+                Self::Big(x) => $ret::Big(x.$method($($arg),*)),
             }
         }
     };
@@ -94,6 +99,7 @@ macro_rules! dispatch_vector_inner {
                 Self::_3(x) => $ret::_3(x.$method($($arg),*)),
                 Self::_5(x) => $ret::_5(x.$method($($arg),*)),
                 Self::_7(x) => $ret::_7(x.$method($($arg),*)),
+                Self::Big(x) => $ret::Big(x.$method($($arg),*)),
             }
         }
     };
@@ -106,6 +112,7 @@ macro_rules! dispatch_vector_inner {
                 Self::_3(x) => $ret::_3(x.$method($($arg),*)),
                 Self::_5(x) => $ret::_5(x.$method($($arg),*)),
                 Self::_7(x) => $ret::_7(x.$method($($arg),*)),
+                Self::Big(x) => $ret::Big(x.$method($($arg),*)),
             }
         }
     };
@@ -118,6 +125,7 @@ macro_rules! dispatch_vector_inner {
                 Self::_3(x) => x.$method($($arg),*),
                 Self::_5(x) => x.$method($($arg),*),
                 Self::_7(x) => x.$method($($arg),*),
+                Self::Big(x) => x.$method($($arg),*),
             }
         }
     };
@@ -129,6 +137,7 @@ macro_rules! dispatch_vector_inner {
                 Self::_3(x) => x.$method($($arg),*),
                 Self::_5(x) => x.$method($($arg),*),
                 Self::_7(x) => x.$method($($arg),*),
+                Self::Big(x) => x.$method($($arg),*),
             }
         }
     };
@@ -140,6 +149,7 @@ macro_rules! dispatch_vector_inner {
                 Self::_3(x) => x.$method($($arg),*),
                 Self::_5(x) => x.$method($($arg),*),
                 Self::_7(x) => x.$method($($arg),*),
+                Self::Big(x) => x.$method($($arg),*),
             }
         }
     }
@@ -155,60 +165,70 @@ macro_rules! dispatch_vector {
     }
 }
 
-macro_rules! match_p {
-    ($p:ident, $($val:tt)*) => {
-        match *$p {
-            2 => Self::_2($($val)*),
-            3 => Self::_3($($val)*),
-            5 => Self::_5($($val)*),
-            7 => Self::_7($($val)*),
-            _ => panic!("Prime not supported: {}", *$p)
-        }
-    }
-}
-
 #[derive(Debug, Hash, Eq, PartialEq, Clone)]
 pub enum FpVector {
-    _2(FpVectorP<2>),
-    _3(FpVectorP<3>),
-    _5(FpVectorP<5>),
-    _7(FpVectorP<7>),
+    _2(FpVectorP<P2>),
+    _3(FpVectorP<P3>),
+    _5(FpVectorP<P5>),
+    _7(FpVectorP<P7>),
+    Big(FpVectorP<ValidPrime>),
 }
 
 #[derive(Debug, Copy, Clone)]
 pub enum Slice<'a> {
-    _2(SliceP<'a, 2>),
-    _3(SliceP<'a, 3>),
-    _5(SliceP<'a, 5>),
-    _7(SliceP<'a, 7>),
+    _2(SliceP<'a, P2>),
+    _3(SliceP<'a, P3>),
+    _5(SliceP<'a, P5>),
+    _7(SliceP<'a, P7>),
+    Big(SliceP<'a, ValidPrime>),
 }
 
 #[derive(Debug)]
 pub enum SliceMut<'a> {
-    _2(SliceMutP<'a, 2>),
-    _3(SliceMutP<'a, 3>),
-    _5(SliceMutP<'a, 5>),
-    _7(SliceMutP<'a, 7>),
+    _2(SliceMutP<'a, P2>),
+    _3(SliceMutP<'a, P3>),
+    _5(SliceMutP<'a, P5>),
+    _7(SliceMutP<'a, P7>),
+    Big(SliceMutP<'a, ValidPrime>),
 }
 
 pub enum FpVectorNonZeroIterator<'a> {
-    _2(FpVectorNonZeroIteratorP<'a, 2>),
-    _3(FpVectorNonZeroIteratorP<'a, 3>),
-    _5(FpVectorNonZeroIteratorP<'a, 5>),
-    _7(FpVectorNonZeroIteratorP<'a, 7>),
+    _2(FpVectorNonZeroIteratorP<'a, P2>),
+    _3(FpVectorNonZeroIteratorP<'a, P3>),
+    _5(FpVectorNonZeroIteratorP<'a, P5>),
+    _7(FpVectorNonZeroIteratorP<'a, P7>),
+    Big(FpVectorNonZeroIteratorP<'a, ValidPrime>),
 }
 
 impl FpVector {
-    pub fn new(p: ValidPrime, len: usize) -> FpVector {
-        match_p!(p, FpVectorP::new_(len))
+    pub fn new<P: Prime>(p: P, len: usize) -> FpVector {
+        match p.as_u32() {
+            2 => FpVector::_2(FpVectorP::new(P2, len)),
+            3 => FpVector::_3(FpVectorP::new(P3, len)),
+            5 => FpVector::_5(FpVectorP::new(P5, len)),
+            7 => FpVector::_7(FpVectorP::new(P7, len)),
+            _ => FpVector::Big(FpVectorP::new(p.to_dyn(), len)),
+        }
     }
 
-    pub fn new_with_capacity(p: ValidPrime, len: usize, capacity: usize) -> FpVector {
-        match_p!(p, FpVectorP::new_with_capacity_(len, capacity))
+    pub fn new_with_capacity<P: Prime>(p: P, len: usize, capacity: usize) -> FpVector {
+        match p.as_u32() {
+            2 => FpVector::_2(FpVectorP::new_with_capacity(P2, len, capacity)),
+            3 => FpVector::_3(FpVectorP::new_with_capacity(P3, len, capacity)),
+            5 => FpVector::_5(FpVectorP::new_with_capacity(P5, len, capacity)),
+            7 => FpVector::_7(FpVectorP::new_with_capacity(P7, len, capacity)),
+            _ => FpVector::Big(FpVectorP::new_with_capacity(p.to_dyn(), len, capacity)),
+        }
     }
 
-    pub fn from_slice(p: ValidPrime, slice: &[u32]) -> Self {
-        match_p!(p, FpVectorP::from(&slice))
+    pub fn from_slice<P: Prime>(p: P, slice: &[u32]) -> Self {
+        match p.as_u32() {
+            2 => FpVector::_2(FpVectorP::from((P2, &slice))),
+            3 => FpVector::_3(FpVectorP::from((P3, &slice))),
+            5 => FpVector::_5(FpVectorP::from((P5, &slice))),
+            7 => FpVector::_7(FpVectorP::from((P7, &slice))),
+            _ => FpVector::Big(FpVectorP::from((p.to_dyn(), &slice))),
+        }
     }
 
     pub fn num_limbs(p: ValidPrime, len: usize) -> usize {
@@ -339,6 +359,7 @@ impl<'a> SliceMut<'a> {
             (SliceMut::_3(x), Slice::_3(y), Slice::_3(z)) => x.add_tensor(offset, coeff, y, z),
             (SliceMut::_5(x), Slice::_5(y), Slice::_5(z)) => x.add_tensor(offset, coeff, y, z),
             (SliceMut::_7(x), Slice::_7(y), Slice::_7(z)) => x.add_tensor(offset, coeff, y, z),
+            (SliceMut::Big(x), Slice::Big(y), Slice::Big(z)) => x.add_tensor(offset, coeff, y, z),
             _ => {
                 panic!("Applying add_tensor to vectors over different primes");
             }
@@ -370,6 +391,7 @@ impl<'a> std::fmt::Display for Slice<'a> {
     fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result {
         if f.alternate() {
             for v in self.iter() {
+                // If self.p >= 11, this will look funky
                 write!(f, "{v}")?;
             }
             Ok(())
@@ -409,7 +431,7 @@ impl<'a> IntoIterator for &'a FpVector {
 }
 
 macro_rules! impl_try_into {
-    ($var:tt, $p:literal) => {
+    ($var:tt, $p:ty) => {
         impl<'a> TryInto<&'a mut FpVectorP<$p>> for &'a mut FpVector {
             type Error = ();
 
@@ -423,10 +445,11 @@ macro_rules! impl_try_into {
     };
 }
 
-impl_try_into!(_2, 2);
-impl_try_into!(_3, 3);
-impl_try_into!(_5, 5);
-impl_try_into!(_7, 7);
+impl_try_into!(_2, P2);
+impl_try_into!(_3, P3);
+impl_try_into!(_5, P5);
+impl_try_into!(_7, P7);
+impl_try_into!(Big, ValidPrime);
 
 impl Serialize for FpVector {
     fn serialize<S>(&self, serializer: S) -> Result<S::Ok, S::Error>
diff --git a/ext/crates/sseq/src/sseq.rs b/ext/crates/sseq/src/sseq.rs
index 53d99174c..ae5823f4c 100644
--- a/ext/crates/sseq/src/sseq.rs
+++ b/ext/crates/sseq/src/sseq.rs
@@ -423,7 +423,7 @@ impl<P: SseqProfile> Sseq<P> {
             }
         }
 
-        let neg_1 = *self.p - 1;
+        let neg_1 = self.p - 1;
 
         let target_r = target_product
             .map(|prod| P::differential_length(x + prod.x - source_x, y + prod.y - source_y))
diff --git a/ext/examples/define_module.rs b/ext/examples/define_module.rs
index 35a8f27c0..80ec8ccd1 100644
--- a/ext/examples/define_module.rs
+++ b/ext/examples/define_module.rs
@@ -9,7 +9,7 @@ use algebra::module::FDModule;
 use algebra::steenrod_evaluator::SteenrodEvaluator;
 use algebra::{AdemAlgebra, Algebra, GeneratedAlgebra};
 use bivec::BiVec;
-use fp::prime::ValidPrime;
+use fp::prime::{Prime, ValidPrime};
 use fp::vector::FpVector;
 
 use anyhow::anyhow;
@@ -86,7 +86,7 @@ pub fn interactive_module_define_fdmodule(
     output_json: &mut Value,
     p: ValidPrime,
 ) -> anyhow::Result<()> {
-    output_json["p"] = Value::from(*p);
+    output_json["p"] = Value::from(p.as_u32());
     let algebra = Arc::new(AdemAlgebra::new(p, false));
 
     let gens = get_gens()?;
@@ -186,7 +186,7 @@ pub fn interactive_module_define_fpmodule(
     }
 
     eprintln!("Input relations");
-    match *p {
+    match p.as_u32() {
         2 => eprintln!("Write relations in the form 'Sq6 * Sq2 * x + Sq7 * y'"),
         _ => eprintln!("Write relations in the form 'Q5 * P(5) * x + 2 * P(1, 3) * Q2 * y', where P(...) and Qi are Milnor basis elements."),
     }
@@ -260,7 +260,7 @@ pub fn interactive_module_define_fpmodule(
         }
     }
 
-    output_json["p"] = Value::from(*p);
+    output_json["p"] = Value::from(p.as_u32());
     output_json["type"] = Value::String("finitely presented module".to_owned());
     for (i, deg_i_gens) in gens.iter_enum() {
         for gen in deg_i_gens {
diff --git a/ext/examples/filtration_one.rs b/ext/examples/filtration_one.rs
index 54594b05b..2f1a0b8e9 100644
--- a/ext/examples/filtration_one.rs
+++ b/ext/examples/filtration_one.rs
@@ -9,7 +9,7 @@ use sseq::coordinates::{Bidegree, BidegreeGenerator};
 
 fn main() -> anyhow::Result<()> {
     let resolution = query_module(None, false)?;
-    assert_eq!(*resolution.prime(), 2);
+    assert_eq!(resolution.prime(), 2);
 
     for b in resolution.iter_stem() {
         let mut i = 0;
diff --git a/ext/examples/massey.rs b/ext/examples/massey.rs
index 8fe665cb1..ab995cbbf 100644
--- a/ext/examples/massey.rs
+++ b/ext/examples/massey.rs
@@ -133,7 +133,7 @@ fn main() -> anyhow::Result<()> {
                 if i != 0 {
                     print!(", ");
                 }
-                print!("{}", entry % *p);
+                print!("{}", entry % p);
             }
             println!("]");
         }
diff --git a/ext/examples/save_bruner.rs b/ext/examples/save_bruner.rs
index 722b83abe..ff216ba3b 100644
--- a/ext/examples/save_bruner.rs
+++ b/ext/examples/save_bruner.rs
@@ -12,7 +12,7 @@ use std::io::{BufWriter, Write as _};
 fn main() -> anyhow::Result<()> {
     let resolution = query_module(Some(AlgebraType::Milnor), false)?;
 
-    assert_eq!(*resolution.prime(), 2);
+    assert_eq!(resolution.prime(), 2);
     let algebra = resolution.algebra();
     let algebra: &MilnorAlgebra = algebra.as_ref().try_into()?;
 
diff --git a/ext/examples/secondary_massey.rs b/ext/examples/secondary_massey.rs
index 1470028a7..981793cdd 100644
--- a/ext/examples/secondary_massey.rs
+++ b/ext/examples/secondary_massey.rs
@@ -438,7 +438,7 @@ fn main() -> anyhow::Result<()> {
             }
             // Now do the -1 part of the null-homotopy of bc.
             {
-                let sign = *p * *p - 1;
+                let sign = p * p - 1;
                 let out = b.product_nullhomotopy(b_tau.as_deref(), &unit_sseq, c, gen);
                 for (i, v) in out.iter_nonzero() {
                     scratch0
@@ -449,11 +449,11 @@ fn main() -> anyhow::Result<()> {
             }
 
             for (i, v) in scratch0.iter().enumerate() {
-                let extra = *v / *p;
-                scratch1.add(&mp[i], extra % *p);
+                let extra = *v / p;
+                scratch1.add(&mp[i], extra % p);
             }
 
-            print!("[{}]", scratch0.iter().map(|x| x % *p).format(", "));
+            print!("[{}]", scratch0.iter().map(|x| *x % p).format(", "));
 
             // Then deal with the rest of the null-homotopy of bc. This is just the null-homotopy
             // of 2.
@@ -467,7 +467,7 @@ fn main() -> anyhow::Result<()> {
                     .for_each(|(a, b)| *a += v * b);
             }
             for (i, v) in scratch0.iter().enumerate() {
-                let extra = (v / *p) % *p;
+                let extra = (*v / p) % p;
                 if extra == 0 {
                     continue;
                 }
diff --git a/ext/examples/secondary_product.rs b/ext/examples/secondary_product.rs
index ca79795c6..d2489bfa1 100644
--- a/ext/examples/secondary_product.rs
+++ b/ext/examples/secondary_product.rs
@@ -83,7 +83,7 @@ fn main() -> anyhow::Result<()> {
             sum.iter_mut().zip_eq(d2).for_each(|(a, b)| *a += x * b);
         }
         assert!(
-            sum.iter().all(|x| x % *p == 0),
+            sum.iter().all(|x| *x % p == 0),
             "Class supports a non-zero d2"
         );
     }
diff --git a/ext/examples/sq0.rs b/ext/examples/sq0.rs
index 10164c4a5..f861f0263 100644
--- a/ext/examples/sq0.rs
+++ b/ext/examples/sq0.rs
@@ -14,7 +14,7 @@ use sseq::coordinates::{Bidegree, BidegreeGenerator};
 fn main() -> anyhow::Result<()> {
     let res = utils::query_module(None, true)?;
     assert!(
-        *res.prime() == 2 && res.target().max_s() == 1 && res.target().module(0).is_unit(),
+        res.prime() == 2 && res.target().max_s() == 1 && res.target().module(0).is_unit(),
         "Sq^0 can only be computed for the sphere at the prime 2"
     );
 
diff --git a/ext/examples/steenrod.rs b/ext/examples/steenrod.rs
index 53ed9a681..dcca5533f 100644
--- a/ext/examples/steenrod.rs
+++ b/ext/examples/steenrod.rs
@@ -19,7 +19,7 @@ fn main() -> anyhow::Result<()> {
     let module = resolution.target().module(0);
     let p = resolution.prime();
 
-    if resolution.target().max_s() != 1 || !module.is_unit() || *p != 2 {
+    if resolution.target().max_s() != 1 || !module.is_unit() || p != 2 {
         panic!("Can only run Steenrod on the sphere");
     }
 
diff --git a/ext/examples/tensor.rs b/ext/examples/tensor.rs
index d97290ff0..1e53135ec 100644
--- a/ext/examples/tensor.rs
+++ b/ext/examples/tensor.rs
@@ -1,7 +1,7 @@
 use algebra::module::{steenrod_module, FDModule, TensorModule};
 use algebra::{AdemAlgebra, SteenrodAlgebra};
 use ext::utils::parse_module_name;
-use fp::prime::ValidPrime;
+use fp::prime::{Prime, ValidPrime};
 
 use anyhow::anyhow;
 use serde_json::json;
@@ -32,7 +32,7 @@ fn main() -> anyhow::Result<()> {
     tensor.name = String::new();
 
     let mut output = json!({
-        "p": *p,
+        "p": p.as_u32(),
     });
     tensor.to_json(&mut output);
 
diff --git a/ext/src/chain_complex/chain_homotopy.rs b/ext/src/chain_complex/chain_homotopy.rs
index bb5da8ec3..0e54e20cd 100644
--- a/ext/src/chain_complex/chain_homotopy.rs
+++ b/ext/src/chain_complex/chain_homotopy.rs
@@ -216,7 +216,7 @@ impl<
 
             self.homotopies[source.s() as i32 - 1].apply(
                 scratch.as_slice_mut(),
-                *p - 1,
+                p - 1,
                 source.t(),
                 self.left
                     .source
diff --git a/ext/src/secondary.rs b/ext/src/secondary.rs
index f538ba9da..19eac270b 100644
--- a/ext/src/secondary.rs
+++ b/ext/src/secondary.rs
@@ -561,7 +561,7 @@ pub trait SecondaryLift: Sync + Sized {
             // Check that we indeed had a lift
             let d = target.differential(target_b.s());
             for (src, tgt) in std::iter::zip(&results, &mut intermediates) {
-                d.apply(tgt.as_slice_mut(), *p - 1, target_b.t(), src.as_slice());
+                d.apply(tgt.as_slice_mut(), p - 1, target_b.t(), src.as_slice());
                 assert!(
                     tgt.is_zero(),
                     "secondary: Failed to lift at {b}. This likely indicates an invalid input."
@@ -864,7 +864,7 @@ where
     }
 
     fn composite(&self, s: u32) -> CompositeData<Self::Algebra> {
-        let p = *self.prime();
+        let p = self.prime();
         // This is -1 mod p^2
         let neg_1 = p * p - 1;
 
@@ -882,7 +882,7 @@ where
 
     fn compute_intermediate(&self, gen: BidegreeGenerator) -> FpVector {
         let p = self.prime();
-        let neg_1 = *p - 1;
+        let neg_1 = p - 1;
         let shifted_b = gen.degree() - self.shift();
         let target = self.target().module(shifted_b.s() - 1);
 
@@ -996,11 +996,11 @@ where
         );
 
         let sign = if (self.underlying.shift.s() as i32 * b.t()) % 2 == 1 {
-            *p * *p - 1
+            p * p - 1
         } else {
             1
         };
-        let filtration_one_sign = if (b.t() % 2) == 1 { *p - 1 } else { 1 };
+        let filtration_one_sign = if (b.t() % 2) == 1 { p - 1 } else { 1 };
 
         let page_data = sseq.map(|sseq| {
             let d = sseq.page_data(tau_source.n(), tau_source.s() as i32);
@@ -1017,14 +1017,14 @@ where
                     .zip_eq(&m0[i])
                     .for_each(|(a, b)| *a += v * b * sign);
                 out.slice_mut(source_num_gens, source_num_gens + tau_num_gens)
-                    .add(m1[i].as_slice(), (v * sign) % *p);
+                    .add(m1[i].as_slice(), (v * sign) % p);
             }
             for (i, v) in scratch0.iter().enumerate() {
-                out.add_basis_element(i, *v % *p);
+                out.add_basis_element(i, *v % p);
 
-                let extra = *v / *p;
+                let extra = *v / p;
                 out.slice_mut(source_num_gens, source_num_gens + tau_num_gens)
-                    .add(mp[i].as_slice(), (extra * filtration_one_sign) % *p);
+                    .add(mp[i].as_slice(), (extra * filtration_one_sign) % p);
             }
             if let Some(page_data) = page_data {
                 page_data.reduce_by_quotient(
@@ -1209,7 +1209,7 @@ where
 
     fn compute_intermediate(&self, gen: BidegreeGenerator) -> FpVector {
         let p = self.prime();
-        let neg_1 = *p - 1;
+        let neg_1 = p - 1;
         let shifted_b = gen.degree() - self.shift();
 
         let target = self.target().module(shifted_b.s() - 1);
@@ -1301,7 +1301,7 @@ where
     }
 
     fn composite(&self, s: u32) -> CompositeData<S::Algebra> {
-        let p = *self.prime();
+        let p = self.prime();
         // This is -1 mod p^2
         let neg_1 = p * p - 1;
 
diff --git a/ext/src/yoneda.rs b/ext/src/yoneda.rs
index bcf072da6..caaa4150c 100644
--- a/ext/src/yoneda.rs
+++ b/ext/src/yoneda.rs
@@ -65,7 +65,7 @@ fn rate_adem_operation(algebra: &AdemAlgebra, deg: i32, idx: usize) -> i32 {
 
 #[allow(dead_code)]
 fn operation_drop(algebra: &AdemAlgebra, deg: i32, idx: usize) -> i32 {
-    if *algebra.prime() != 2 {
+    if algebra.prime() != 2 {
         return 1;
     }
     let elt = algebra.basis_element_from_index(deg, idx);
diff --git a/web_ext/sseq_gui/src/managers.rs b/web_ext/sseq_gui/src/managers.rs
index 62da7cb69..cd7783953 100644
--- a/web_ext/sseq_gui/src/managers.rs
+++ b/web_ext/sseq_gui/src/managers.rs
@@ -127,7 +127,7 @@ impl ResolutionManager {
             let unit_resolution = Resolution::new_from_json(
                 json!({
                     "type": "finite dimensional module",
-                    "p": *resolution.prime(),
+                    "p": resolution.prime(),
                     "gens": {"x0": 0},
                     "actions": [],
                 }),
diff --git a/web_ext/sseq_gui/src/resolution_wrapper.rs b/web_ext/sseq_gui/src/resolution_wrapper.rs
index 79e149f00..aa4ec1050 100644
--- a/web_ext/sseq_gui/src/resolution_wrapper.rs
+++ b/web_ext/sseq_gui/src/resolution_wrapper.rs
@@ -190,7 +190,7 @@ impl<CC: ChainComplex> Resolution<CC> {
             || (!left && mult.t() * (source.s() as i32) % 2 != 0)
         {
             for entry in product.iter_mut().flatten() {
-                *entry = ((*p - 1) * *entry) % *p;
+                *entry = ((p - 1) * *entry) % p;
             }
         }
 
@@ -317,7 +317,7 @@ impl<CC: ChainComplex> Resolution<CC> {
                         val += elt.class[i] * result.entry(idx);
                     }
                 }
-                products[k].push(val % *self.prime());
+                products[k].push(val % self.prime());
             }
         }
         self.add_structline(&elt.name, source, elt.b, true, products);


Test patch:
diff --git a/ext/crates/fp/tests/row_reduce.rs b/ext/crates/fp/tests/row_reduce.rs
index 37a008124..a8f539a37 100644
--- a/ext/crates/fp/tests/row_reduce.rs
+++ b/ext/crates/fp/tests/row_reduce.rs
@@ -1,15 +1,36 @@
+use std::sync::OnceLock;
+
 use proptest::prelude::*;
 
-use fp::{matrix::Matrix, prime::ValidPrime};
+use fp::{
+    matrix::Matrix,
+    prime::{Prime, ValidPrime},
+};
 
-/// An arbitrary `ValidPrime`
+/// An arbitrary `ValidPrime` in the range `2..(1 << 24)`, plus the largest prime that we support.
 fn arb_prime() -> impl Strategy<Value = ValidPrime> {
-    prop_oneof![
-        Just(ValidPrime::new(2)),
-        Just(ValidPrime::new(3)),
-        Just(ValidPrime::new(5)),
-        Just(ValidPrime::new(7)),
-    ]
+    static TEST_PRIMES: OnceLock<Vec<ValidPrime>> = OnceLock::new();
+    let test_primes = TEST_PRIMES.get_or_init(|| {
+        // Sieve of erathosthenes
+        const MAX: usize = 1 << 24;
+        let mut is_prime = Vec::new();
+        is_prime.resize_with(MAX, || true);
+        is_prime[0] = false;
+        is_prime[1] = false;
+        for i in 2..MAX {
+            if is_prime[i] {
+                for j in ((2 * i)..MAX).step_by(i) {
+                    is_prime[j] = false;
+                }
+            }
+        }
+        (0..MAX)
+            .filter(|&i| is_prime[i])
+            .map(|p| ValidPrime::new_unchecked(p as u32))
+            .chain(std::iter::once(ValidPrime::new_unchecked(2147483647)))
+            .collect()
+    });
+    (0..test_primes.len()).prop_map(|i| test_primes[i])
 }
 
 /// An increasing sequence of numbers between 0 and `cols`, where the sequence has a length between
@@ -46,7 +67,7 @@ fn arb_coeff_row_pair_seq(
     p: ValidPrime,
     rows: usize,
 ) -> impl Strategy<Value = Vec<(u32, (usize, usize))>> {
-    proptest::collection::vec((1..*p, arb_row_pair(rows)), 0..1000)
+    proptest::collection::vec((1..p.as_u32(), arb_row_pair(rows)), 0..1000)
 }
 
 // This is a macro used to define functions that take in values from strategies as arguments.