diff --git a/src/g1.rs b/src/g1.rs
index b211761..d72b192 100644
--- a/src/g1.rs
+++ b/src/g1.rs
@@ -7,6 +7,7 @@ use core::{
     borrow::Borrow,
     fmt,
     iter::Sum,
+    mem::transmute,
     ops::{Add, AddAssign, Mul, MulAssign, Neg, Sub, SubAssign},
 };
 use elliptic_curve::consts::U48;
@@ -530,6 +531,16 @@ impl G1Affine {
             ))
         }
     }
+
+    pub fn sum_of_products(points: &[Self], scalars: &[Scalar]) -> G1Projective {
+        let num_components = points.len().min(scalars.len());
+        let mut raw_scalars = Vec::with_capacity(num_components * 32);
+        for s in scalars {
+            raw_scalars.extend_from_slice(&s.to_le_bytes());
+        }
+        let inner_points = unsafe { transmute::<_, &[blst_p1_affine]>(points) };
+        G1Projective(inner_points.mult(&raw_scalars, 255))
+    }
 }
 
 /// This is an element of $\mathbb{G}_1$ represented in the projective coordinate space.
@@ -1642,6 +1653,29 @@ mod tests {
         assert_eq!(naive, pippenger);
     }
 
+    #[test]
+    fn test_multi_exp_affine() {
+        const SIZE: usize = 10;
+        let mut rng = XorShiftRng::from_seed([
+            0x59, 0x62, 0xbe, 0x5d, 0x76, 0x3d, 0x31, 0x8d, 0x17, 0xdb, 0x37, 0x32, 0x54, 0x06,
+            0xbc, 0xe5,
+        ]);
+
+        let points: Vec<G1Affine> = (0..SIZE)
+            .map(|_| G1Projective::random(&mut rng).to_affine())
+            .collect();
+        let scalars: Vec<Scalar> = (0..SIZE).map(|_| Scalar::random(&mut rng)).collect();
+
+        let mut naive = points[0] * scalars[0];
+        for i in 1..SIZE {
+            naive += points[i] * scalars[i];
+        }
+
+        let pippenger = G1Affine::sum_of_products(points.as_slice(), scalars.as_slice());
+
+        assert_eq!(naive, pippenger);
+    }
+
     #[test]
     fn test_identity() {
         let id = G1Projective::IDENTITY;
diff --git a/src/g2.rs b/src/g2.rs
index f63a677..4791f0c 100644
--- a/src/g2.rs
+++ b/src/g2.rs
@@ -7,6 +7,7 @@ use core::{
     borrow::Borrow,
     fmt,
     iter::Sum,
+    mem::transmute,
     ops::{Add, AddAssign, Mul, MulAssign, Neg, Sub, SubAssign},
 };
 use std::hash::Hash;
@@ -506,6 +507,16 @@ impl G2Affine {
     pub const fn compressed_size() -> usize {
         COMPRESSED_SIZE
     }
+
+    pub fn sum_of_products(points: &[Self], scalars: &[Scalar]) -> G2Projective {
+        let num_components = points.len().min(scalars.len());
+        let mut raw_scalars = Vec::with_capacity(num_components * 32);
+        for s in scalars {
+            raw_scalars.extend_from_slice(&s.to_le_bytes());
+        }
+        let inner_points = unsafe { transmute::<_, &[blst_p2_affine]>(points) };
+        G2Projective(inner_points.mult(&raw_scalars, 255))
+    }
 }
 
 /// This is an element of $\mathbb{G}_2$ represented in the projective coordinate space.
@@ -1603,6 +1614,29 @@ mod tests {
         assert_eq!(naive, pippenger);
     }
 
+    #[test]
+    fn test_multi_exp_affine() {
+        const SIZE: usize = 10;
+        let mut rng = XorShiftRng::from_seed([
+            0x59, 0x62, 0xbe, 0x5d, 0x76, 0x3d, 0x31, 0x8d, 0x17, 0xdb, 0x37, 0x32, 0x54, 0x06,
+            0xbc, 0xe5,
+        ]);
+
+        let points: Vec<G2Affine> = (0..SIZE)
+            .map(|_| G2Projective::random(&mut rng).to_affine())
+            .collect();
+        let scalars: Vec<Scalar> = (0..SIZE).map(|_| Scalar::random(&mut rng)).collect();
+
+        let mut naive = points[0] * scalars[0];
+        for i in 1..SIZE {
+            naive += points[i] * scalars[i];
+        }
+
+        let pippenger = G2Affine::sum_of_products(points.as_slice(), scalars.as_slice());
+
+        assert_eq!(naive, pippenger);
+    }
+
     #[test]
     fn test_hex() {
         let g1 = G2Projective::GENERATOR;