diff --git a/src/c_kzg_4844.c b/src/c_kzg_4844.c
index 65e28d2..b843d43 100644
--- a/src/c_kzg_4844.c
+++ b/src/c_kzg_4844.c
@@ -92,7 +92,7 @@ static const g2_t g2_generator = {{{{
  * The first 32 roots of unity in the finite field F_r.
  *
  * For element `{A, B, C, D}`, the field element value is `A + B * 2^64 + C * 2^128 + D * 2^192`. This format may be
- * converted to an `fr_t` type via the #fr_from_uint64s library function.
+ * converted to an `fr_t` type via the #blst_fr_from_uint64 function.
  *
  * The decimal values may be calculated with the following Python code:
  * @code{.py}
@@ -235,16 +235,6 @@ static int log_2_byte(byte b) {
     return r;
 }
 
-/**
- * Create a field element from an array of four 64-bit unsigned integers.
- *
- * @param out  The field element equivalent of @p vals
- * @param vals The array of 64-bit integers to be converted, little-endian ordering of the 64-bit words
- */
-static void fr_from_uint64s(fr_t *out, const uint64_t vals[4]) {
-    blst_fr_from_uint64(out, vals);
-}
-
 /**
  * Test whether the operand is one in the finite field.
  *
@@ -275,39 +265,6 @@ static bool fr_equal(const fr_t *aa, const fr_t *bb) {
     return a[0] == b[0] && a[1] == b[1] && a[2] == b[2] && a[3] == b[3];
 }
 
-/**
- * Add two field elements.
- *
- * @param[out] out @p a plus @p b in the field
- * @param[in]  a   Field element
- * @param[in]  b   Field element
- */
-static void fr_add(fr_t *out, const fr_t *a, const fr_t *b) {
-    blst_fr_add(out, a, b);
-}
-
-/**
- * Subtract one field element from another.
- *
- * @param[out] out @p a minus @p b in the field
- * @param[in]  a   Field element
- * @param[in]  b   Field element
- */
-static void fr_sub(fr_t *out, const fr_t *a, const fr_t *b) {
-    blst_fr_sub(out, a, b);
-}
-
-/**
- * Multiply two field elements.
- *
- * @param[out] out @p a multiplied by @p b in the field
- * @param[in]  a   Multiplicand
- * @param[in]  b   Multiplier
- */
-static void fr_mul(fr_t *out, const fr_t *a, const fr_t *b) {
-    blst_fr_mul(out, a, b);
-}
-
 /**
  * Divide a field element by another.
  *
@@ -321,16 +278,6 @@ static void fr_div(fr_t *out, const fr_t *a, const fr_t *b) {
     blst_fr_mul(out, a, &tmp);
 }
 
-/**
- * Square a field element.
- *
- * @param[out] out @p a squared
- * @param[in]  a   A field element
- */
-static void fr_sqr(fr_t *out, const fr_t *a) {
-    blst_fr_sqr(out, a);
-}
-
 /**
  * Exponentiation of a field element.
  *
@@ -348,10 +295,10 @@ static void fr_pow(fr_t *out, const fr_t *a, uint64_t n) {
 
     while (true) {
         if (n & 1) {
-            fr_mul(out, out, &tmp);
+            blst_fr_mul(out, out, &tmp);
         }
         if ((n >>= 1) == 0) break;
-        fr_sqr(&tmp, &tmp);
+        blst_fr_sqr(&tmp, &tmp);
     }
 }
 
@@ -365,17 +312,7 @@ static void fr_pow(fr_t *out, const fr_t *a, uint64_t n) {
  */
 static void fr_from_uint64(fr_t *out, uint64_t n) {
     uint64_t vals[] = {n, 0, 0, 0};
-    fr_from_uint64s(out, vals);
-}
-
-/**
- * Inverse of a field element.
- *
- * @param[out] out The inverse of @p a
- * @param[in]  a   A field element
- */
-static void fr_inv(fr_t *out, const fr_t *a) {
-    blst_fr_eucl_inverse(out, a);
+    blst_fr_from_uint64(out, vals);
 }
 
 /**
@@ -397,14 +334,14 @@ static C_KZG_RET fr_batch_inv(fr_t *out, const fr_t *a, size_t len) {
     prod[0] = a[0];
 
     for(i = 1; i < len; i++) {
-        fr_mul(&prod[i], &a[i], &prod[i - 1]);
+        blst_fr_mul(&prod[i], &a[i], &prod[i - 1]);
     }
 
     blst_fr_eucl_inverse(&inv, &prod[len - 1]);
 
     for(i = len - 1; i > 0; i--) {
-        fr_mul(&out[i], &inv, &prod[i - 1]);
-        fr_mul(&inv, &a[i], &inv);
+        blst_fr_mul(&out[i], &inv, &prod[i - 1]);
+        blst_fr_mul(&inv, &a[i], &inv);
     }
     out[0] = inv;
 
@@ -413,19 +350,6 @@ out:
     return ret;
 }
 
-/**
- * Add or double G1 points.
- *
- * This is safe if the two points are the same.
- *
- * @param[out] out @p a plus @p b in the group
- * @param[in]  a   G1 group point
- * @param[in]  b   G1 group point
- */
-static void g1_add_or_dbl(g1_t *out, const g1_t *a, const g1_t *b) {
-    blst_p1_add_or_double(out, a, b);
-}
-
 /**
  * Multiply a G1 group element by a field element.
  *
@@ -837,8 +761,8 @@ static void poly_lincomb(Polynomial *out, const Polynomial *vectors, const fr_t
         out->evals[j] = fr_zero;
     for (i = 0; i < n; i++) {
         for (j = 0; j < FIELD_ELEMENTS_PER_BLOB; j++) {
-            fr_mul(&tmp, &scalars[i], &vectors[i].evals[j]);
-            fr_add(&out->evals[j], &out->evals[j], &tmp);
+            blst_fr_mul(&tmp, &scalars[i], &vectors[i].evals[j]);
+            blst_fr_add(&out->evals[j], &out->evals[j], &tmp);
         }
     }
 }
@@ -848,7 +772,7 @@ static void compute_powers(fr_t *out, fr_t *x, uint64_t n) {
     fr_t current_power = fr_one;
     for (uint64_t i = 0; i < n; i++) {
         out[i] = current_power;
-        fr_mul(&current_power, &current_power, x);
+        blst_fr_mul(&current_power, &current_power, x);
     }
 }
 
@@ -875,7 +799,7 @@ static C_KZG_RET evaluate_polynomial_in_evaluation_form(fr_t *out, const Polynom
             ret = C_KZG_OK;
             goto out;
         }
-        fr_sub(&inverses_in[i], x, &roots_of_unity[i]);
+        blst_fr_sub(&inverses_in[i], x, &roots_of_unity[i]);
     }
 
     ret = fr_batch_inv(inverses, inverses_in, FIELD_ELEMENTS_PER_BLOB);
@@ -883,15 +807,15 @@ static C_KZG_RET evaluate_polynomial_in_evaluation_form(fr_t *out, const Polynom
 
     *out = fr_zero;
     for (i = 0; i < FIELD_ELEMENTS_PER_BLOB; i++) {
-        fr_mul(&tmp, &inverses[i], &roots_of_unity[i]);
-        fr_mul(&tmp, &tmp, &p->evals[i]);
-        fr_add(out, out, &tmp);
+        blst_fr_mul(&tmp, &inverses[i], &roots_of_unity[i]);
+        blst_fr_mul(&tmp, &tmp, &p->evals[i]);
+        blst_fr_add(out, out, &tmp);
     }
     fr_from_uint64(&tmp, FIELD_ELEMENTS_PER_BLOB);
     fr_div(out, out, &tmp);
     fr_pow(&tmp, x, FIELD_ELEMENTS_PER_BLOB);
-    fr_sub(&tmp, &tmp, &fr_one);
-    fr_mul(out, out, &tmp);
+    blst_fr_sub(&tmp, &tmp, &fr_one);
+    blst_fr_mul(out, out, &tmp);
 
 out:
     free(inverses_in);
@@ -1008,15 +932,15 @@ static C_KZG_RET compute_kzg_proof(g1_t *out, const Polynomial *p, const fr_t *x
             continue;
         }
         // (p_i - y) / (ω_i - x)
-        fr_sub(&q.evals[i], &p->evals[i], &y);
-        fr_sub(&inverses_in[i], &roots_of_unity[i], x);
+        blst_fr_sub(&q.evals[i], &p->evals[i], &y);
+        blst_fr_sub(&inverses_in[i], &roots_of_unity[i], x);
     }
 
     ret = fr_batch_inv(inverses, inverses_in, FIELD_ELEMENTS_PER_BLOB);
     if (ret != C_KZG_OK) goto out;
 
     for (i = 0; i < FIELD_ELEMENTS_PER_BLOB; i++) {
-        fr_mul(&q.evals[i], &q.evals[i], &inverses[i]);
+        blst_fr_mul(&q.evals[i], &q.evals[i], &inverses[i]);
     }
 
     if (m) { // ω_m == x
@@ -1024,16 +948,16 @@ static C_KZG_RET compute_kzg_proof(g1_t *out, const Polynomial *p, const fr_t *x
         for (i = 0; i < FIELD_ELEMENTS_PER_BLOB; i++) {
             if (i == m) continue;
             // (p_i - y) * ω_i / (x * (x - ω_i))
-            fr_sub(&tmp, x, &roots_of_unity[i]);
-            fr_mul(&inverses_in[i], &tmp, x);
+            blst_fr_sub(&tmp, x, &roots_of_unity[i]);
+            blst_fr_mul(&inverses_in[i], &tmp, x);
         }
         ret = fr_batch_inv(inverses, inverses_in, FIELD_ELEMENTS_PER_BLOB);
         if (ret != C_KZG_OK) goto out;
         for (i = 0; i < FIELD_ELEMENTS_PER_BLOB; i++) {
-            fr_sub(&tmp, &p->evals[i], &y);
-            fr_mul(&tmp, &tmp, &roots_of_unity[i]);
-            fr_mul(&tmp, &tmp, &inverses[i]);
-            fr_add(&q.evals[m], &q.evals[m], &tmp);
+            blst_fr_sub(&tmp, &p->evals[i], &y);
+            blst_fr_mul(&tmp, &tmp, &roots_of_unity[i]);
+            blst_fr_mul(&tmp, &tmp, &inverses[i]);
+            blst_fr_add(&q.evals[m], &q.evals[m], &tmp);
         }
     }
 
@@ -1195,7 +1119,7 @@ static void fft_g1_fast(g1_t *out, const g1_t *in, uint64_t stride, const fr_t *
             g1_t y_times_root;
             g1_mul(&y_times_root, &out[i + half], &roots[i * roots_stride]);
             g1_sub(&out[i + half], &out[i], &y_times_root);
-            g1_add_or_dbl(&out[i], &out[i], &y_times_root);
+            blst_p1_add_or_double(&out[i], &out[i], &y_times_root);
         }
     } else {
         *out = *in;
@@ -1220,7 +1144,7 @@ static C_KZG_RET fft_g1(g1_t *out, const g1_t *in, bool inverse, uint64_t n, con
     if (inverse) {
         fr_t inv_len;
         fr_from_uint64(&inv_len, n);
-        fr_inv(&inv_len, &inv_len);
+        blst_fr_eucl_inverse(&inv_len, &inv_len);
         fft_g1_fast(out, in, 1, fs->reverse_roots_of_unity, stride, n);
         for (uint64_t i = 0; i < n; i++) {
             g1_mul(&out[i], &out[i], &inv_len);
@@ -1249,7 +1173,7 @@ static C_KZG_RET expand_root_of_unity(fr_t *out, const fr_t *root, uint64_t widt
 
     for (uint64_t i = 2; !fr_is_one(&out[i - 1]); i++) {
         CHECK(i <= width);
-        fr_mul(&out[i], &out[i - 1], root);
+        blst_fr_mul(&out[i], &out[i - 1], root);
     }
     CHECK(fr_is_one(&out[width]));
 
@@ -1286,7 +1210,7 @@ static C_KZG_RET new_fft_settings(FFTSettings *fs, unsigned int max_scale) {
     fs->roots_of_unity = NULL;
 
     CHECK((max_scale < sizeof scale2_root_of_unity / sizeof scale2_root_of_unity[0]));
-    fr_from_uint64s(&root_of_unity, scale2_root_of_unity[max_scale]);
+    blst_fr_from_uint64(&root_of_unity, scale2_root_of_unity[max_scale]);
 
     // Allocate space for the roots of unity
     ret = new_fr_array(&fs->expanded_roots_of_unity, fs->max_width + 1);