diff --git a/crypto/src/commitments/kzg.rs b/crypto/src/commitments/kzg.rs
index 2ecbb483c..c65d64d8c 100644
--- a/crypto/src/commitments/kzg.rs
+++ b/crypto/src/commitments/kzg.rs
@@ -1,11 +1,16 @@
 use super::traits::IsCommitmentScheme;
 use alloc::{borrow::ToOwned, vec::Vec};
 use core::{marker::PhantomData, mem};
+use lambdaworks_math::gpu::icicle::GpuMSMPoint;
+use lambdaworks_math::traits::ByteConversion;
 use lambdaworks_math::{
     cyclic_group::IsGroup,
     elliptic_curve::traits::IsPairing,
     errors::DeserializationError,
-    field::{element::FieldElement, traits::IsPrimeField},
+    field::{
+        element::FieldElement,
+        traits::{IsField, IsPrimeField},
+    },
     msm::pippenger::msm,
     polynomial::Polynomial,
     traits::{AsBytes, Deserializable},
@@ -136,12 +141,12 @@ where
 }
 
 #[derive(Clone)]
-pub struct KateZaveruchaGoldberg<F: IsPrimeField, P: IsPairing> {
+pub struct KateZaveruchaGoldberg<F: IsField, P: IsPairing> {
     srs: StructuredReferenceString<P::G1Point, P::G2Point>,
     phantom: PhantomData<F>,
 }
 
-impl<F: IsPrimeField, P: IsPairing> KateZaveruchaGoldberg<F, P> {
+impl<F: IsField, P: IsPairing> KateZaveruchaGoldberg<F, P> {
     pub fn new(srs: StructuredReferenceString<P::G1Point, P::G2Point>) -> Self {
         Self {
             srs,
@@ -150,20 +155,22 @@ impl<F: IsPrimeField, P: IsPairing> KateZaveruchaGoldberg<F, P> {
     }
 }
 
-impl<const N: usize, F: IsPrimeField<RepresentativeType = UnsignedInteger<N>>, P: IsPairing>
-    IsCommitmentScheme<F> for KateZaveruchaGoldberg<F, P>
+impl<
+        const N: usize,
+        F: IsField<BaseType = UnsignedInteger<N>>
+            + IsPrimeField<RepresentativeType = UnsignedInteger<N>>,
+        P: IsPairing,
+    > IsCommitmentScheme<F> for KateZaveruchaGoldberg<F, P>
+where
+    FieldElement<F>: ByteConversion,
+    P::G1Point: GpuMSMPoint,
 {
     type Commitment = P::G1Point;
 
     fn commit(&self, p: &Polynomial<FieldElement<F>>) -> Self::Commitment {
-        let coefficients: Vec<_> = p
-            .coefficients
-            .iter()
-            .map(|coefficient| coefficient.representative())
-            .collect();
         msm(
-            &coefficients,
-            &self.srs.powers_main_group[..coefficients.len()],
+            &p.coefficients,
+            &self.srs.powers_main_group[..p.coefficients.len()],
         )
         .expect("`points` is sliced by `cs`'s length")
     }
diff --git a/math/Cargo.toml b/math/Cargo.toml
index 02df065d4..847afb37b 100644
--- a/math/Cargo.toml
+++ b/math/Cargo.toml
@@ -24,11 +24,17 @@ objc = { version = "0.2.7", optional = true }
 # cuda
 cudarc = { version = "0.9.7", optional = true }
 
+# Icicle integration
+icicle-cuda-runtime = { git = "https://github.com/ingonyama-zk/icicle.git", tag = "v1.4.0", optional = true }
+icicle-core = { git = "https://github.com/ingonyama-zk/icicle.git", tag = "v1.4.0", optional = true }
+icicle-bls12-377 = { git = "https://github.com/ingonyama-zk/icicle.git", tag = "v1.4.0", optional = true}
+icicle-bls12-381 = { git = "https://github.com/ingonyama-zk/icicle.git", tag = "v1.4.0", optional = true }
+icicle-bn254 = { git = "https://github.com/ingonyama-zk/icicle.git", tag = "v1.4.0", optional = true }
 
 lambdaworks-gpu = { workspace = true, optional = true }
 
 [dev-dependencies]
-rand = { version = "0.8.5", default-features = false }
+rand = { version = "0.8.5" }
 rand_chacha = "0.3.1"
 criterion = "0.5.1"
 const-random = "0.1.15"
@@ -45,6 +51,7 @@ lambdaworks-serde-binary = ["dep:serde", "alloc"]
 lambdaworks-serde-string = ["dep:serde", "dep:serde_json", "alloc"]
 proptest = ["dep:proptest"]
 winter_compatibility = ["winter-math", "miden-core"]
+icicle = ["dep:icicle-cuda-runtime", "dep:icicle-core", "dep:icicle-bls12-377", "dep:icicle-bls12-381", "dep:icicle-bn254"]
 
 # gpu
 metal = [
@@ -84,6 +91,11 @@ name = "criterion_msm"
 harness = false
 required-features = ["parallel"]
 
+[[bench]]
+name = "criterion_icicle"
+harness = false
+required-features = ["icicle"]
+
 [[bench]]
 name = "criterion_fft"
 harness = false
diff --git a/math/benches/criterion_icicle.rs b/math/benches/criterion_icicle.rs
new file mode 100644
index 000000000..07a271220
--- /dev/null
+++ b/math/benches/criterion_icicle.rs
@@ -0,0 +1,71 @@
+use criterion::{black_box, criterion_group, criterion_main, Criterion};
+use lambdaworks_math::{
+    cyclic_group::IsGroup,
+    elliptic_curve::{
+        short_weierstrass::curves::{
+            bls12_377::curve::BLS12377Curve, bls12_381::curve::BLS12381Curve,
+            bn_254::curve::BN254Curve,
+        },
+        traits::IsEllipticCurve,
+    },
+    field::{element::FieldElement, traits::IsField},
+};
+
+use lambdaworks_math::gpu::icicle::{
+    bls12_377::bls12_377_g1_msm, bls12_381::bls12_381_g1_msm, bn254::bn254_g1_msm,
+};
+use rand::{rngs::StdRng, Rng, SeedableRng};
+
+pub fn generate_cs_and_points<C: IsEllipticCurve>(
+    msm_size: usize,
+) -> (Vec<FieldElement<C::BaseField>>, Vec<C::PointRepresentation>)
+where
+    <C::BaseField as IsField>::BaseType: From<u64>,
+{
+    // We use a seeded rng so the benchmarks are reproducible.
+    let mut rng = StdRng::seed_from_u64(42);
+
+    let g = C::generator();
+
+    let cs: Vec<_> = (0..msm_size)
+        .map(|_| FieldElement::<C::BaseField>::new(rng.gen::<u64>().into()))
+        .collect();
+
+    let points: Vec<_> = (0..msm_size)
+        .map(|_| g.operate_with_self(rng.gen::<u64>()))
+        .collect();
+
+    (cs, points)
+}
+
+pub fn msm_benchmarks_with_size(c: &mut Criterion, msm_size: usize) {
+    let mut group = c.benchmark_group(format!("MSM benchmarks with size {msm_size}"));
+
+    let (cs, points) = generate_cs_and_points::<BLS12381Curve>(msm_size);
+    group.bench_function("BLS12_381", |bench| {
+        bench.iter(|| black_box(bls12_381_g1_msm(&cs, &points, None)));
+    });
+
+    let (cs, points) = generate_cs_and_points::<BLS12377Curve>(msm_size);
+    group.bench_function("BLS12_377", |bench| {
+        bench.iter(|| black_box(bls12_377_g1_msm(&cs, &points, None)));
+    });
+
+    let (cs, points) = generate_cs_and_points::<BN254Curve>(msm_size);
+    group.bench_function("BN_254", |bench| {
+        bench.iter(|| black_box(bn254_g1_msm(&cs, &points, None)));
+    });
+}
+
+pub fn run_benchmarks(c: &mut Criterion) {
+    let exponents = 1..=18;
+
+    for exp in exponents {
+        let msm_size = 1 << exp;
+
+        msm_benchmarks_with_size(c, msm_size);
+    }
+}
+
+criterion_group!(icicle_msm, run_benchmarks);
+criterion_main!(icicle_msm);
diff --git a/math/src/elliptic_curve/short_weierstrass/curves/bls12_377/curve.rs b/math/src/elliptic_curve/short_weierstrass/curves/bls12_377/curve.rs
index c324b35de..a9b2958b1 100644
--- a/math/src/elliptic_curve/short_weierstrass/curves/bls12_377/curve.rs
+++ b/math/src/elliptic_curve/short_weierstrass/curves/bls12_377/curve.rs
@@ -5,6 +5,8 @@ use crate::{
     elliptic_curve::short_weierstrass::traits::IsShortWeierstrass, field::element::FieldElement,
 };
 
+pub type BLS12377FieldElement = FieldElement<BLS12377PrimeField>;
+
 /// The description of the curve.
 #[derive(Clone, Debug)]
 pub struct BLS12377Curve;
diff --git a/math/src/fft/errors.rs b/math/src/fft/errors.rs
index 4bbec91c2..e27e7d74d 100644
--- a/math/src/fft/errors.rs
+++ b/math/src/fft/errors.rs
@@ -8,6 +8,9 @@ use lambdaworks_gpu::metal::abstractions::errors::MetalError;
 #[cfg(feature = "cuda")]
 use lambdaworks_gpu::cuda::abstractions::errors::CudaError;
 
+#[cfg(feature = "icicle")]
+use icicle_core::error::IcicleError;
+
 #[derive(Debug)]
 pub enum FFTError {
     RootOfUnityError(u64),
@@ -17,6 +20,8 @@ pub enum FFTError {
     MetalError(MetalError),
     #[cfg(feature = "cuda")]
     CudaError(CudaError),
+    #[cfg(feature = "icicle")]
+    IcicleError(IcicleError),
 }
 
 impl Display for FFTError {
@@ -37,6 +42,8 @@ impl Display for FFTError {
             FFTError::CudaError(_) => {
                 write!(f, "A CUDA related error has ocurred")
             }
+            #[cfg(feature = "icicle")]
+            FFTError::IcicleError(e) => write!(f, "Icicle GPU backend failure. {:?}", e),
         }
     }
 }
diff --git a/math/src/fft/polynomial.rs b/math/src/fft/polynomial.rs
index a27811d31..cbb341a3c 100644
--- a/math/src/fft/polynomial.rs
+++ b/math/src/fft/polynomial.rs
@@ -7,6 +7,8 @@ use crate::{
         traits::{IsFFTField, RootsConfig},
     },
     polynomial::Polynomial,
+    traits::ByteConversion,
+    gpu::icicle::{IcicleFFT, GpuMSMPoint}
 };
 use alloc::{vec, vec::Vec};
 
@@ -14,10 +16,15 @@ use alloc::{vec, vec::Vec};
 use crate::fft::gpu::cuda::polynomial::{evaluate_fft_cuda, interpolate_fft_cuda};
 #[cfg(feature = "metal")]
 use crate::fft::gpu::metal::polynomial::{evaluate_fft_metal, interpolate_fft_metal};
+#[cfg(feature = "icicle")]
+use crate::gpu::icicle::{evaluate_fft_icicle, interpolate_fft_icicle};
 
 use super::cpu::{ops, roots_of_unity};
 
-impl<E: IsField> Polynomial<FieldElement<E>> {
+impl<E: IsField + IsFFTField + IcicleFFT> Polynomial<FieldElement<E>>
+where
+    FieldElement<E>: ByteConversion,
+{
     /// Returns `N` evaluations of this polynomial using FFT over a domain in a subfield F of E (so the results
     /// are P(w^i), with w being a primitive root of unity).
     /// `N = max(self.coeff_len(), domain_size).next_power_of_two() * blowup_factor`.
@@ -26,7 +33,10 @@ impl<E: IsField> Polynomial<FieldElement<E>> {
         poly: &Polynomial<FieldElement<E>>,
         blowup_factor: usize,
         domain_size: Option<usize>,
-    ) -> Result<Vec<FieldElement<E>>, FFTError> {
+    ) -> Result<Vec<FieldElement<E>>, FFTError>
+    where
+        FieldElement<E>: ByteConversion,
+    {
         let domain_size = domain_size.unwrap_or(0);
         let len = core::cmp::max(poly.coeff_len(), domain_size).next_power_of_two() * blowup_factor;
 
@@ -51,17 +61,20 @@ impl<E: IsField> Polynomial<FieldElement<E>> {
             }
         }
 
-        #[cfg(feature = "cuda")]
+        #[cfg(feature = "icicle")]
         {
-            // TODO: support multiple fields with CUDA
-            if F::field_name() == "stark256" {
-                Ok(evaluate_fft_cuda(&coeffs)?)
+            if !F::field_name().is_empty() {
+                Ok(evaluate_fft_icicle::<F, E>(&coeffs)?)
             } else {
+                println!(
+                    "GPU evaluation failed for field {}. Program will fallback to CPU.",
+                    core::any::type_name::<F>()
+                );
                 evaluate_fft_cpu::<F, E>(&coeffs)
             }
         }
 
-        #[cfg(all(not(feature = "metal"), not(feature = "cuda")))]
+        #[cfg(all(not(feature = "metal"), not(feature = "icicle")))]
         {
             evaluate_fft_cpu::<F, E>(&coeffs)
         }
@@ -76,7 +89,10 @@ impl<E: IsField> Polynomial<FieldElement<E>> {
         blowup_factor: usize,
         domain_size: Option<usize>,
         offset: &FieldElement<F>,
-    ) -> Result<Vec<FieldElement<E>>, FFTError> {
+    ) -> Result<Vec<FieldElement<E>>, FFTError>
+    where
+        FieldElement<E>: ByteConversion,
+    {
         let scaled = poly.scale(offset);
         Polynomial::evaluate_fft::<F>(&scaled, blowup_factor, domain_size)
     }
@@ -100,16 +116,20 @@ impl<E: IsField> Polynomial<FieldElement<E>> {
             }
         }
 
-        #[cfg(feature = "cuda")]
+        #[cfg(feature = "icicle")]
         {
             if !F::field_name().is_empty() {
-                Ok(interpolate_fft_cuda(fft_evals)?)
+                Ok(interpolate_fft_icicle::<F, E>(fft_evals)?)
             } else {
-                interpolate_fft_cpu::<F, E>(fft_evals)
+                println!(
+                    "GPU evaluation failed for field {}. Program will fallback to CPU.",
+                    core::any::type_name::<F>()
+                );
+                interpolate_fft_cpu::<F, E>(&fft_evals)
             }
         }
 
-        #[cfg(all(not(feature = "metal"), not(feature = "cuda")))]
+        #[cfg(all(not(feature = "metal"), not(feature = "icicle")))]
         {
             interpolate_fft_cpu::<F, E>(fft_evals)
         }
@@ -127,13 +147,13 @@ impl<E: IsField> Polynomial<FieldElement<E>> {
     }
 }
 
-pub fn compose_fft<F, E>(
-    poly_1: &Polynomial<FieldElement<E>>,
-    poly_2: &Polynomial<FieldElement<E>>,
-) -> Polynomial<FieldElement<E>>
+pub fn compose_fft<F>(
+    poly_1: &Polynomial<FieldElement<F>>,
+    poly_2: &Polynomial<FieldElement<F>>,
+) -> Polynomial<FieldElement<F>>
 where
-    F: IsFFTField + IsSubFieldOf<E>,
-    E: IsField,
+    F: IsFFTField + IcicleFFT,
+    FieldElement<F>: ByteConversion,
 {
     let poly_2_evaluations = Polynomial::evaluate_fft::<F>(poly_2, 1, None).unwrap();
 
diff --git a/math/src/gpu/icicle/bls12_377.rs b/math/src/gpu/icicle/bls12_377.rs
new file mode 100644
index 000000000..2381e11be
--- /dev/null
+++ b/math/src/gpu/icicle/bls12_377.rs
@@ -0,0 +1,156 @@
+use icicle_bls12_377::curve;
+use icicle_core::{error::IcicleError, msm, traits::FieldImpl};
+use icicle_cuda_runtime::{memory::HostOrDeviceSlice, stream::CudaStream};
+
+use crate::{
+    elliptic_curve::short_weierstrass::{
+        curves::bls12_377::{
+            curve::{BLS12377Curve, BLS12377FieldElement},
+            field_extension::BLS12377PrimeField,
+        },
+        point::ShortWeierstrassProjectivePoint,
+    },
+    errors::ByteConversionError,
+    field::element::FieldElement,
+    traits::ByteConversion,
+};
+
+impl BLS12377FieldElement {
+    fn to_icicle(&self) -> curve::BaseField {
+        curve::BaseField::from_bytes_le(&self.to_bytes_le())
+    }
+
+    fn to_icicle_scalar(&self) -> curve::ScalarField {
+        curve::ScalarField::from_bytes_le(&self.to_bytes_le())
+    }
+
+    fn from_icicle(icicle: &curve::BaseField) -> Result<Self, ByteConversionError> {
+        Self::from_bytes_le(&icicle.to_bytes_le())
+    }
+}
+
+impl ShortWeierstrassProjectivePoint<BLS12377Curve> {
+    fn to_icicle(&self) -> curve::G1Affine {
+        let s = self.to_affine();
+        curve::G1Affine {
+            x: s.x().to_icicle(),
+            y: s.y().to_icicle(),
+        }
+    }
+
+    fn from_icicle(icicle: &curve::G1Projective) -> Result<Self, ByteConversionError> {
+        Ok(Self::new([
+            FieldElement::<BLS12377PrimeField>::from_icicle(&icicle.x).unwrap(),
+            FieldElement::<BLS12377PrimeField>::from_icicle(&icicle.y).unwrap(),
+            FieldElement::<BLS12377PrimeField>::from_icicle(&icicle.z).unwrap(),
+        ]))
+    }
+}
+
+pub fn bls12_377_g1_msm(
+    scalars: &[BLS12377FieldElement],
+    points: &[ShortWeierstrassProjectivePoint<BLS12377Curve>],
+    config: Option<msm::MSMConfig>,
+) -> Result<ShortWeierstrassProjectivePoint<BLS12377Curve>, IcicleError> {
+    let mut cfg = config.unwrap_or(msm::MSMConfig::default());
+    let scalars = HostOrDeviceSlice::Host(
+        scalars
+            .iter()
+            .map(|scalar| scalar.to_icicle_scalar())
+            .collect::<Vec<_>>(),
+    );
+    let points = HostOrDeviceSlice::Host(
+        points
+            .iter()
+            .map(|point| point.to_icicle())
+            .collect::<Vec<_>>(),
+    );
+    let mut msm_results = HostOrDeviceSlice::cuda_malloc(1).unwrap();
+    let stream = CudaStream::create().unwrap();
+    cfg.ctx.stream = &stream;
+    cfg.is_async = true;
+    msm::msm(&scalars, &points, &cfg, &mut msm_results).unwrap();
+    let mut msm_host_result = vec![curve::G1Projective::zero(); 1];
+    stream.synchronize().unwrap();
+    msm_results.copy_to_host(&mut msm_host_result[..]).unwrap();
+    stream.destroy().unwrap();
+    let res =
+        ShortWeierstrassProjectivePoint::<BLS12377Curve>::from_icicle(&msm_host_result[0]).unwrap();
+    Ok(res)
+}
+
+#[cfg(test)]
+mod test {
+    use super::*;
+    use crate::{
+        elliptic_curve::{
+            short_weierstrass::curves::bls12_377::curve::BLS12377FieldElement,
+            traits::IsEllipticCurve,
+        },
+        field::element::FieldElement,
+        msm::pippenger::msm,
+    };
+
+    impl ShortWeierstrassProjectivePoint<BLS12377Curve> {
+        fn from_icicle_affine(
+            icicle: &curve::G1Affine,
+        ) -> Result<ShortWeierstrassProjectivePoint<BLS12377Curve>, ByteConversionError> {
+            Ok(Self::new([
+                FieldElement::<BLS12377PrimeField>::from_icicle(&icicle.x).unwrap(),
+                FieldElement::<BLS12377PrimeField>::from_icicle(&icicle.y).unwrap(),
+                FieldElement::one(),
+            ]))
+        }
+    }
+
+    fn point_times_5() -> ShortWeierstrassProjectivePoint<BLS12377Curve> {
+        let x = BLS12377FieldElement::from_hex_unchecked(
+            "3c852d5aab73fbb51e57fbf5a0a8b5d6513ec922b2611b7547bfed74cba0dcdfc3ad2eac2733a4f55d198ec82b9964",
+        );
+        let y = BLS12377FieldElement::from_hex_unchecked(
+            "a71425e68e55299c64d7eada9ae9c3fb87a9626b941d17128b64685fc07d0e635f3c3a512903b4e0a43e464045967b",
+        );
+        BLS12377Curve::create_point_from_affine(x, y).unwrap()
+    }
+
+    #[test]
+    fn to_from_icicle() {
+        // convert value of 5 to icicle and back again and that icicle 5 matches
+        let point = point_times_5();
+        let icicle_point = point.to_icicle();
+        let res =
+            ShortWeierstrassProjectivePoint::<BLS12377Curve>::from_icicle_affine(&icicle_point)
+                .unwrap();
+        assert_eq!(point, res)
+    }
+
+    #[test]
+    fn to_from_icicle_generator() {
+        // Convert generator and see that it matches
+        let point = BLS12377Curve::generator();
+        let icicle_point = point.to_icicle();
+        let res =
+            ShortWeierstrassProjectivePoint::<BLS12377Curve>::from_icicle_affine(&icicle_point)
+                .unwrap();
+        assert_eq!(point, res)
+    }
+
+    #[test]
+    fn icicle_g1_msm() {
+        const LEN: usize = 20;
+        let eight: BLS12377FieldElement = FieldElement::from(8);
+        let lambda_scalars = vec![eight; LEN];
+        let lambda_points = (0..LEN).map(|_| point_times_5()).collect::<Vec<_>>();
+        let expected = msm(
+            &lambda_scalars
+                .clone()
+                .into_iter()
+                .map(|x| x.representative())
+                .collect::<Vec<_>>(),
+            &lambda_points,
+        )
+        .unwrap();
+        let res = bls12_377_g1_msm(&lambda_scalars, &lambda_points, None).unwrap();
+        assert_eq!(res, expected);
+    }
+}
diff --git a/math/src/gpu/icicle/bls12_381.rs b/math/src/gpu/icicle/bls12_381.rs
new file mode 100644
index 000000000..ce1f123b8
--- /dev/null
+++ b/math/src/gpu/icicle/bls12_381.rs
@@ -0,0 +1,152 @@
+use icicle_bls12_381::curve;
+use icicle_core::{error::IcicleError, msm, traits::FieldImpl};
+use icicle_cuda_runtime::{memory::HostOrDeviceSlice, stream::CudaStream};
+
+use crate::{
+    elliptic_curve::short_weierstrass::{
+        curves::bls12_381::{
+            curve::{BLS12381Curve, BLS12381FieldElement},
+            field_extension::BLS12381PrimeField,
+        },
+        point::ShortWeierstrassProjectivePoint,
+    },
+    errors::ByteConversionError,
+    field::element::FieldElement,
+    traits::ByteConversion,
+};
+
+impl BLS12381FieldElement {
+    fn to_icicle(&self) -> curve::BaseField {
+        curve::BaseField::from_bytes_le(&self.to_bytes_le())
+    }
+
+    fn to_icicle_scalar(&self) -> curve::ScalarField {
+        curve::ScalarField::from_bytes_le(&self.to_bytes_le())
+    }
+
+    fn from_icicle(icicle: &curve::BaseField) -> Result<Self, ByteConversionError> {
+        Self::from_bytes_le(&icicle.to_bytes_le())
+    }
+}
+
+impl ShortWeierstrassProjectivePoint<BLS12381Curve> {
+    fn to_icicle(&self) -> curve::G1Affine {
+        let s = self.to_affine();
+        curve::G1Affine {
+            x: s.x().to_icicle(),
+            y: s.y().to_icicle(),
+        }
+    }
+
+    fn from_icicle(icicle: &curve::G1Projective) -> Result<Self, ByteConversionError> {
+        Ok(Self::new([
+            FieldElement::<BLS12381PrimeField>::from_icicle(&icicle.x).unwrap(),
+            FieldElement::<BLS12381PrimeField>::from_icicle(&icicle.y).unwrap(),
+            FieldElement::<BLS12381PrimeField>::from_icicle(&icicle.z).unwrap(),
+        ]))
+    }
+}
+
+pub fn bls12_381_g1_msm(
+    scalars: &[BLS12381FieldElement],
+    points: &[ShortWeierstrassProjectivePoint<BLS12381Curve>],
+    config: Option<msm::MSMConfig>,
+) -> Result<ShortWeierstrassProjectivePoint<BLS12381Curve>, IcicleError> {
+    let mut cfg = config.unwrap_or(msm::MSMConfig::default());
+    let scalars = HostOrDeviceSlice::Host(
+        scalars
+            .iter()
+            .map(|scalar| scalar.to_icicle_scalar())
+            .collect::<Vec<_>>(),
+    );
+    let points = HostOrDeviceSlice::Host(
+        points
+            .iter()
+            .map(|point| point.to_icicle())
+            .collect::<Vec<_>>(),
+    );
+    let mut msm_results = HostOrDeviceSlice::cuda_malloc(1).unwrap();
+    let stream = CudaStream::create().unwrap();
+    cfg.ctx.stream = &stream;
+    cfg.is_async = true;
+    msm::msm(&scalars, &points, &cfg, &mut msm_results).unwrap();
+    let mut msm_host_result = vec![curve::G1Projective::zero(); 1];
+    stream.synchronize().unwrap();
+    msm_results.copy_to_host(&mut msm_host_result[..]).unwrap();
+    stream.destroy().unwrap();
+    let res =
+        ShortWeierstrassProjectivePoint::<BLS12381Curve>::from_icicle(&msm_host_result[0]).unwrap();
+    Ok(res)
+}
+
+#[cfg(test)]
+mod test {
+    use super::*;
+    use crate::{
+        elliptic_curve::{
+            short_weierstrass::curves::bls12_381::curve::BLS12381FieldElement,
+            traits::IsEllipticCurve,
+        },
+        field::element::FieldElement,
+        msm::pippenger::msm,
+    };
+
+    impl ShortWeierstrassProjectivePoint<BLS12381Curve> {
+        fn from_icicle_affine(
+            icicle: &curve::G1Affine,
+        ) -> Result<ShortWeierstrassProjectivePoint<BLS12381Curve>, ByteConversionError> {
+            Ok(Self::new([
+                FieldElement::<BLS12381PrimeField>::from_icicle(&icicle.x).unwrap(),
+                FieldElement::<BLS12381PrimeField>::from_icicle(&icicle.y).unwrap(),
+                FieldElement::one(),
+            ]))
+        }
+    }
+
+    fn point_times_5() -> ShortWeierstrassProjectivePoint<BLS12381Curve> {
+        let x = BLS12381FieldElement::from_hex_unchecked(
+            "32bcce7e71eb50384918e0c9809f73bde357027c6bf15092dd849aa0eac274d43af4c68a65fb2cda381734af5eecd5c",
+        );
+        let y = BLS12381FieldElement::from_hex_unchecked(
+            "11e48467b19458aabe7c8a42dc4b67d7390fdf1e150534caadddc7e6f729d8890b68a5ea6885a21b555186452b954d88",
+        );
+        BLS12381Curve::create_point_from_affine(x, y).unwrap()
+    }
+
+    #[test]
+    fn to_from_icicle() {
+        // convert value of 5 to icicle and back again and that icicle 5 matches
+        let point = point_times_5();
+        let icicle_point = point.to_icicle();
+        let res =
+            ShortWeierstrassProjectivePoint::<BLS12381Curve>::from_icicle_affine(&icicle_point)
+                .unwrap();
+        assert_eq!(point, res)
+    }
+
+    #[test]
+    fn to_from_icicle_generator() {
+        // Convert generator and see that it matches
+        let point = BLS12381Curve::generator();
+        let icicle_point = point.to_icicle();
+        let res =
+            ShortWeierstrassProjectivePoint::<BLS12381Curve>::from_icicle_affine(&icicle_point)
+                .unwrap();
+        assert_eq!(point, res)
+    }
+
+    #[test]
+    fn icicle_g1_msm() {
+        const LEN: usize = 20;
+        let eight: BLS12381FieldElement = FieldElement::from(8);
+        let lambda_scalars = vec![eight; LEN];
+        let lambda_points = (0..LEN).map(|_| point_times_5()).collect::<Vec<_>>();
+        let expected = msm(
+            &lambda_scalars,
+            &lambda_points,
+        )
+        .unwrap();
+        let res = bls12_381_g1_msm(&lambda_scalars, &lambda_points, None).unwrap();
+        assert_eq!(res, expected);
+    }
+}
diff --git a/math/src/gpu/icicle/bn254.rs b/math/src/gpu/icicle/bn254.rs
new file mode 100644
index 000000000..daf91a61f
--- /dev/null
+++ b/math/src/gpu/icicle/bn254.rs
@@ -0,0 +1,153 @@
+use icicle_bn254::curve;
+use icicle_core::{error::IcicleError, msm, traits::FieldImpl};
+use icicle_cuda_runtime::{memory::HostOrDeviceSlice, stream::CudaStream};
+
+use crate::{
+    elliptic_curve::short_weierstrass::{
+        curves::bn_254::{
+            curve::{BN254Curve, BN254FieldElement},
+            field_extension::BN254PrimeField,
+        },
+        point::ShortWeierstrassProjectivePoint,
+    },
+    errors::ByteConversionError,
+    field::element::FieldElement,
+    traits::ByteConversion,
+};
+
+impl BN254FieldElement {
+    fn to_icicle(&self) -> curve::BaseField {
+        curve::BaseField::from_bytes_le(&self.to_bytes_le())
+    }
+
+    fn to_icicle_scalar(&self) -> curve::ScalarField {
+        curve::ScalarField::from_bytes_le(&self.to_bytes_le())
+    }
+
+    fn from_icicle(icicle: &curve::BaseField) -> Result<Self, ByteConversionError> {
+        Self::from_bytes_le(&icicle.to_bytes_le())
+    }
+}
+
+impl ShortWeierstrassProjectivePoint<BN254Curve> {
+    fn to_icicle(&self) -> curve::G1Affine {
+        let s = self.to_affine();
+        curve::G1Affine {
+            x: s.x().to_icicle(),
+            y: s.y().to_icicle(),
+        }
+    }
+
+    fn from_icicle(icicle: &curve::G1Projective) -> Result<Self, ByteConversionError> {
+        Ok(Self::new([
+            FieldElement::<BN254PrimeField>::from_icicle(&icicle.x).unwrap(),
+            FieldElement::<BN254PrimeField>::from_icicle(&icicle.y).unwrap(),
+            FieldElement::<BN254PrimeField>::from_icicle(&icicle.z).unwrap(),
+        ]))
+    }
+}
+
+pub fn bn254_g1_msm(
+    scalars: &[BN254FieldElement],
+    points: &[ShortWeierstrassProjectivePoint<BN254Curve>],
+    config: Option<msm::MSMConfig>,
+) -> Result<ShortWeierstrassProjectivePoint<BN254Curve>, IcicleError> {
+    let mut cfg = config.unwrap_or(msm::MSMConfig::default());
+    let scalars = HostOrDeviceSlice::Host(
+        scalars
+            .iter()
+            .map(|scalar| scalar.to_icicle_scalar())
+            .collect::<Vec<_>>(),
+    );
+    let points = HostOrDeviceSlice::Host(
+        points
+            .iter()
+            .map(|point| point.to_icicle())
+            .collect::<Vec<_>>(),
+    );
+    let mut msm_results = HostOrDeviceSlice::cuda_malloc(1).unwrap();
+    let stream = CudaStream::create().unwrap();
+    cfg.ctx.stream = &stream;
+    cfg.is_async = true;
+    msm::msm(&scalars, &points, &cfg, &mut msm_results).unwrap();
+    let mut msm_host_result = vec![curve::G1Projective::zero(); 1];
+    stream.synchronize().unwrap();
+    msm_results.copy_to_host(&mut msm_host_result[..]).unwrap();
+    stream.destroy().unwrap();
+    let res =
+        ShortWeierstrassProjectivePoint::<BN254Curve>::from_icicle(&msm_host_result[0]).unwrap();
+    Ok(res)
+}
+
+#[cfg(test)]
+mod test {
+    use super::*;
+    use crate::{
+        elliptic_curve::{
+            short_weierstrass::curves::bn_254::curve::BN254FieldElement, traits::IsEllipticCurve,
+        },
+        field::element::FieldElement,
+        msm::pippenger::msm,
+    };
+
+    impl ShortWeierstrassProjectivePoint<BN254Curve> {
+        fn from_icicle_affine(
+            icicle: &curve::G1Affine,
+        ) -> Result<ShortWeierstrassProjectivePoint<BN254Curve>, ByteConversionError> {
+            Ok(Self::new([
+                FieldElement::<BN254PrimeField>::from_icicle(&icicle.x).unwrap(),
+                FieldElement::<BN254PrimeField>::from_icicle(&icicle.y).unwrap(),
+                FieldElement::one(),
+            ]))
+        }
+    }
+
+    fn point_times_5() -> ShortWeierstrassProjectivePoint<BN254Curve> {
+        let x = BN254FieldElement::from_hex_unchecked(
+            "16ab03b69dfb4f870b0143ebf6a71b7b2e4053ca7a4421d09a913b8b834bbfa3",
+        );
+        let y = BN254FieldElement::from_hex_unchecked(
+            "2512347279ba1049ef97d4ec348d838f939d2b7623e88f4826643cf3889599b2",
+        );
+        BN254Curve::create_point_from_affine(x, y).unwrap()
+    }
+
+    #[test]
+    fn to_from_icicle() {
+        // convert value of 5 to icicle and back again and that icicle 5 matches
+        let point = point_times_5();
+        let icicle_point = point.to_icicle();
+        let res = ShortWeierstrassProjectivePoint::<BN254Curve>::from_icicle_affine(&icicle_point)
+            .unwrap();
+        assert_eq!(point, res)
+    }
+
+    #[test]
+    fn to_from_icicle_generator() {
+        // Convert generator and see that it matches
+        let point = BN254Curve::generator();
+        let icicle_point = point.to_icicle();
+        let res = ShortWeierstrassProjectivePoint::<BN254Curve>::from_icicle_affine(&icicle_point)
+            .unwrap();
+        assert_eq!(point, res)
+    }
+
+    #[test]
+    fn icicle_g1_msm() {
+        const LEN: usize = 20;
+        let eight: BN254FieldElement = FieldElement::from(8);
+        let lambda_scalars = vec![eight; LEN];
+        let lambda_points = (0..LEN).map(|_| point_times_5()).collect::<Vec<_>>();
+        let expected = msm(
+            &lambda_scalars
+                .clone()
+                .into_iter()
+                .map(|x| x.representative())
+                .collect::<Vec<_>>(),
+            &lambda_points,
+        )
+        .unwrap();
+        let res = bn254_g1_msm(&lambda_scalars, &lambda_points, None).unwrap();
+        assert_eq!(res, expected);
+    }
+}
diff --git a/math/src/gpu/icicle/mod.rs b/math/src/gpu/icicle/mod.rs
new file mode 100644
index 000000000..170bce439
--- /dev/null
+++ b/math/src/gpu/icicle/mod.rs
@@ -0,0 +1,398 @@
+use crate::{
+    cyclic_group::IsGroup,
+    elliptic_curve::{
+        short_weierstrass::{
+            curves::{
+                bls12_377::curve::BLS12377Curve,
+                bls12_381::{
+                    curve::BLS12381Curve, default_types::FrField,
+                    field_extension::BLS12381PrimeField, twist::BLS12381TwistCurve,
+                },
+                bn_254::{curve::BN254Curve, field_extension::BN254PrimeField},
+            },
+            point::ShortWeierstrassProjectivePoint,
+        },
+        traits::IsEllipticCurve,
+    },
+    errors::ByteConversionError,
+    fft::errors::FFTError,
+    field::{
+        element::FieldElement,
+        fields::fft_friendly::stark_252_prime_field::Stark252PrimeField,
+        traits::{IsFFTField, IsField, IsSubFieldOf},
+    },
+    msm::naive::MSMError,
+    polynomial::Polynomial,
+    traits::ByteConversion,
+};
+use icicle_bls12_377::curve::CurveCfg as IcicleBLS12377Curve;
+use icicle_bls12_381::curve::{
+    CurveCfg as IcicleBLS12381Curve, ScalarCfg as IcicleBLS12381ScalarCfg,
+};
+use icicle_bn254::curve::{CurveCfg as IcicleBN254Curve, ScalarCfg as IcicleBN254ScalarCfg};
+use icicle_core::{
+    curve::{Affine, Curve, Projective},
+    msm,
+    ntt::{NTTConfig, NTTDir, NTT},
+    traits::FieldImpl,
+};
+use icicle_cuda_runtime::{memory::HostOrDeviceSlice, stream::CudaStream};
+
+use std::fmt::Debug;
+
+impl GpuMSMPoint for ShortWeierstrassProjectivePoint<BLS12381Curve> {
+    type LambdaCurve = BLS12381Curve;
+    type GpuCurve = IcicleBLS12381Curve;
+
+    fn curve_name() -> &'static str {
+        "BLS12381"
+    }
+
+    fn to_icicle_affine(point: &Self) -> Affine<Self::GpuCurve> {
+        let s = Self::to_affine(point);
+        Affine::<Self::GpuCurve> {
+            x: Self::to_icicle_field(s.x()),
+            y: Self::to_icicle_field(s.y()),
+        }
+    }
+
+    fn from_icicle_projective(
+        icicle: &Projective<Self::GpuCurve>,
+    ) -> Result<Self, ByteConversionError> {
+        Ok(Self::new([
+            Self::from_icicle_field(&icicle.x).unwrap(),
+            Self::from_icicle_field(&icicle.y).unwrap(),
+            Self::from_icicle_field(&icicle.z).unwrap(),
+        ]))
+    }
+}
+
+//NOTE THIS IS A PLACEHOLDER COMPILING ICICLE G2 TOOK TO LONG!
+impl GpuMSMPoint for ShortWeierstrassProjectivePoint<BLS12381TwistCurve> {
+    type LambdaCurve = BLS12381Curve;
+    type GpuCurve = IcicleBLS12381Curve;
+
+    fn curve_name() -> &'static str {
+        ""
+    }
+
+    fn to_icicle_affine(point: &Self) -> Affine<Self::GpuCurve> {
+        let s = Self::to_affine(point);
+        Affine::<Self::GpuCurve> {
+            x: Self::to_icicle_field(s.x()),
+            y: Self::to_icicle_field(s.y()),
+        }
+    }
+
+    fn from_icicle_projective(
+        icicle: &Projective<Self::GpuCurve>,
+    ) -> Result<Self, ByteConversionError> {
+        Ok(Self::new([
+            Self::from_icicle_field(&icicle.x).unwrap(),
+            Self::from_icicle_field(&icicle.y).unwrap(),
+            Self::from_icicle_field(&icicle.z).unwrap(),
+        ]))
+    }
+}
+
+impl GpuMSMPoint for ShortWeierstrassProjectivePoint<BLS12377Curve> {
+    type LambdaCurve = BLS12377Curve;
+    type GpuCurve = IcicleBLS12377Curve;
+    fn curve_name() -> &'static str {
+        "BLS12377"
+    }
+
+    fn to_icicle_affine(point: &Self) -> Affine<Self::GpuCurve> {
+        let s = Self::to_affine(point);
+        Affine::<Self::GpuCurve> {
+            x: Self::to_icicle_field(s.x()),
+            y: Self::to_icicle_field(s.y()),
+        }
+    }
+
+    fn from_icicle_projective(
+        icicle: &Projective<Self::GpuCurve>,
+    ) -> Result<Self, ByteConversionError> {
+        Ok(Self::new([
+            Self::from_icicle_field(&icicle.x).unwrap(),
+            Self::from_icicle_field(&icicle.y).unwrap(),
+            Self::from_icicle_field(&icicle.z).unwrap(),
+        ]))
+    }
+}
+
+impl GpuMSMPoint for ShortWeierstrassProjectivePoint<BN254Curve> {
+    type LambdaCurve = BN254Curve;
+    type GpuCurve = IcicleBN254Curve;
+    fn curve_name() -> &'static str {
+        "BN254"
+    }
+
+    fn to_icicle_affine(point: &Self) -> Affine<Self::GpuCurve> {
+        let s = Self::to_affine(point);
+        Affine::<Self::GpuCurve> {
+            x: Self::to_icicle_field(s.x()),
+            y: Self::to_icicle_field(s.y()),
+        }
+    }
+
+    fn from_icicle_projective(
+        icicle: &Projective<Self::GpuCurve>,
+    ) -> Result<Self, ByteConversionError> {
+        Ok(Self::new([
+            Self::from_icicle_field(&icicle.x).unwrap(),
+            Self::from_icicle_field(&icicle.y).unwrap(),
+            Self::from_icicle_field(&icicle.z).unwrap(),
+        ]))
+    }
+}
+
+pub trait GpuMSMPoint: IsGroup {
+    type LambdaCurve: IsEllipticCurve + Clone + Debug;
+    type GpuCurve: Curve + msm::MSM<Self::GpuCurve>;
+    //type FE: ByteConversion;
+    /// Used for searching this field's implementation in other languages, e.g in MSL
+    /// for executing parallel operations with the Metal API.
+    fn curve_name() -> &'static str {
+        ""
+    }
+
+    fn to_icicle_affine(point: &Self) -> Affine<Self::GpuCurve>;
+
+    fn from_icicle_projective(
+        icicle: &Projective<Self::GpuCurve>,
+    ) -> Result<Self, ByteConversionError>;
+
+    fn to_icicle_field<FE: ByteConversion>(element: &FE) -> <Self::GpuCurve as Curve>::BaseField {
+        <Self::GpuCurve as Curve>::BaseField::from_bytes_le(&element.to_bytes_le())
+    }
+
+    fn to_icicle_scalar<FE: ByteConversion>(
+        element: &FE,
+    ) -> <Self::GpuCurve as Curve>::ScalarField {
+        <Self::GpuCurve as Curve>::ScalarField::from_bytes_le(&element.to_bytes_le())
+    }
+
+    fn from_icicle_field<FE: ByteConversion>(
+        icicle: &<Self::GpuCurve as Curve>::BaseField,
+    ) -> Result<FE, ByteConversionError> {
+        FE::from_bytes_le(&icicle.to_bytes_le())
+    }
+}
+
+pub trait IcicleFFT: IsField
+where
+    FieldElement<Self>: ByteConversion,
+{
+    type ScalarField: FieldImpl;
+    type Config: NTT<<Self as IcicleFFT>::ScalarField>;
+
+    fn to_icicle_scalar(element: &FieldElement<Self>) -> Self::ScalarField {
+        Self::ScalarField::from_bytes_le(&element.to_bytes_le())
+    }
+
+    fn from_icicle_scalar(
+        icicle: &Self::ScalarField,
+    ) -> Result<FieldElement<Self>, ByteConversionError> {
+        FieldElement::<Self>::from_bytes_le(&icicle.to_bytes_le())
+    }
+}
+
+impl IcicleFFT for BLS12381PrimeField {
+    type ScalarField = <IcicleBLS12381Curve as Curve>::ScalarField;
+    type Config = IcicleBLS12381ScalarCfg;
+}
+
+impl IcicleFFT for FrField {
+    type ScalarField = <IcicleBLS12381Curve as Curve>::ScalarField;
+    type Config = IcicleBLS12381ScalarCfg;
+}
+
+// DUMMY IMPLEMENTATION OF STARK252 -> Fails when Icicle feature flag enabled
+impl IcicleFFT for Stark252PrimeField {
+    type ScalarField = <IcicleBLS12381Curve as Curve>::ScalarField;
+    type Config = IcicleBLS12381ScalarCfg;
+}
+
+impl IcicleFFT for BN254PrimeField {
+    type ScalarField = <IcicleBN254Curve as Curve>::ScalarField;
+    type Config = IcicleBN254ScalarCfg;
+}
+
+pub fn icicle_msm<F: IsField, G: GpuMSMPoint>(
+    cs: &[FieldElement<F>],
+    points: &[G],
+) -> Result<G, MSMError>
+where
+    FieldElement<F>: ByteConversion,
+{
+    let mut cfg = msm::MSMConfig::default();
+    let scalars = HostOrDeviceSlice::Host(
+        cs.iter()
+            .map(|scalar| G::to_icicle_scalar(scalar))
+            .collect::<Vec<_>>(),
+    );
+    let points = HostOrDeviceSlice::Host(
+        points
+            .iter()
+            .map(|point| G::to_icicle_affine(point))
+            .collect::<Vec<_>>(),
+    );
+    let mut msm_results = HostOrDeviceSlice::cuda_malloc(1).unwrap();
+    let stream = CudaStream::create().unwrap();
+    cfg.ctx.stream = &stream;
+    cfg.is_async = true;
+    msm::msm(&scalars, &points, &cfg, &mut msm_results).unwrap();
+    let mut msm_host_result = [Projective::<G::GpuCurve>::zero(); 1];
+    stream.synchronize().unwrap();
+    msm_results.copy_to_host(&mut msm_host_result[..]).unwrap();
+    stream.destroy().unwrap();
+    let res = G::from_icicle_projective(&msm_host_result[0]).unwrap();
+    Ok(res)
+}
+
+pub fn evaluate_fft_icicle<F, E>(
+    coeffs: &Vec<FieldElement<E>>,
+) -> Result<Vec<FieldElement<E>>, FFTError>
+where
+    F: IsFFTField + IsSubFieldOf<E>,
+    FieldElement<E>: ByteConversion,
+    E: IsField + IcicleFFT,
+{
+    let size = coeffs.len();
+    let mut cfg = NTTConfig::default();
+    let order = coeffs.len() as u64;
+    let dir = NTTDir::kForward;
+    let scalars = HostOrDeviceSlice::Host(
+        coeffs
+            .iter()
+            .map(|scalar| E::to_icicle_scalar(&scalar))
+            .collect::<Vec<_>>(),
+    );
+    let mut ntt_results = HostOrDeviceSlice::cuda_malloc(size).unwrap();
+    let stream = CudaStream::create().unwrap();
+    cfg.ctx.stream = &stream;
+    cfg.is_async = true;
+    let root_of_unity = E::to_icicle_scalar(
+        &(F::get_primitive_root_of_unity(order).unwrap() * FieldElement::<E>::one()),
+    );
+    <E as IcicleFFT>::Config::initialize_domain(root_of_unity, &cfg.ctx).unwrap();
+    <E as IcicleFFT>::Config::ntt_unchecked(&scalars, dir, &cfg, &mut ntt_results).unwrap();
+    stream.synchronize().unwrap();
+    let mut ntt_host_results = vec![E::ScalarField::zero(); size];
+    ntt_results.copy_to_host(&mut ntt_host_results[..]).unwrap();
+    stream.destroy().unwrap();
+    let res = ntt_host_results
+        .as_slice()
+        .iter()
+        .map(|scalar| E::from_icicle_scalar(&scalar).unwrap())
+        .collect::<Vec<_>>();
+    Ok(res)
+}
+
+pub fn interpolate_fft_icicle<F, E>(
+    coeffs: &[FieldElement<E>],
+) -> Result<Polynomial<FieldElement<E>>, FFTError>
+where
+    F: IsFFTField + IsSubFieldOf<E>,
+    FieldElement<E>: ByteConversion,
+    E: IsField + IcicleFFT,
+{
+    let size = coeffs.len();
+    let mut cfg = NTTConfig::default();
+    let order = coeffs.len() as u64;
+    let dir = NTTDir::kInverse;
+    let scalars = HostOrDeviceSlice::Host(
+        coeffs
+            .iter()
+            .map(|scalar| E::to_icicle_scalar(scalar))
+            .collect::<Vec<_>>(),
+    );
+    let mut ntt_results = HostOrDeviceSlice::cuda_malloc(size).unwrap();
+    let stream = CudaStream::create().unwrap();
+    cfg.ctx.stream = &stream;
+    cfg.is_async = true;
+    let root_of_unity = E::to_icicle_scalar(
+        &(F::get_primitive_root_of_unity(order).unwrap() * FieldElement::<E>::one()),
+    );
+    <E as IcicleFFT>::Config::initialize_domain(root_of_unity, &cfg.ctx).unwrap();
+    <E as IcicleFFT>::Config::ntt_unchecked(&scalars, dir, &cfg, &mut ntt_results).unwrap();
+    stream.synchronize().unwrap();
+    let mut ntt_host_results = vec![E::ScalarField::zero(); size];
+    ntt_results.copy_to_host(&mut ntt_host_results[..]).unwrap();
+    stream.destroy().unwrap();
+    let res = ntt_host_results
+        .as_slice()
+        .iter()
+        .map(|scalar| E::from_icicle_scalar(&scalar).unwrap())
+        .collect::<Vec<_>>();
+    Ok(Polynomial::new(&res))
+}
+
+#[cfg(test)]
+mod test {
+    use super::*;
+    use crate::{
+        elliptic_curve::{
+            short_weierstrass::curves::bls12_381::curve::BLS12381FieldElement,
+            traits::IsEllipticCurve,
+        },
+        field::element::FieldElement,
+        msm::pippenger::msm,
+    };
+
+    impl ShortWeierstrassProjectivePoint<BLS12381Curve> {
+        fn from_icicle_affine(
+            icicle: &curve::G1Affine,
+        ) -> Result<ShortWeierstrassProjectivePoint<BLS12381Curve>, ByteConversionError> {
+            Ok(Self::new([
+                FieldElement::<BLS12381PrimeField>::from_icicle(&icicle.x).unwrap(),
+                FieldElement::<BLS12381PrimeField>::from_icicle(&icicle.y).unwrap(),
+                FieldElement::one(),
+            ]))
+        }
+    }
+
+    fn point_times_5() -> ShortWeierstrassProjectivePoint<BLS12381Curve> {
+        let x = BLS12381FieldElement::from_hex_unchecked(
+            "32bcce7e71eb50384918e0c9809f73bde357027c6bf15092dd849aa0eac274d43af4c68a65fb2cda381734af5eecd5c",
+        );
+        let y = BLS12381FieldElement::from_hex_unchecked(
+            "11e48467b19458aabe7c8a42dc4b67d7390fdf1e150534caadddc7e6f729d8890b68a5ea6885a21b555186452b954d88",
+        );
+        BLS12381Curve::create_point_from_affine(x, y).unwrap()
+    }
+
+    #[test]
+    fn to_from_icicle() {
+        // convert value of 5 to icicle and back again and that icicle 5 matches
+        let point = point_times_5();
+        let icicle_point = point.to_icicle();
+        let res =
+            ShortWeierstrassProjectivePoint::<BLS12381Curve>::from_icicle_affine(&icicle_point)
+                .unwrap();
+        assert_eq!(point, res)
+    }
+
+    #[test]
+    fn to_from_icicle_generator() {
+        // Convert generator and see that it matches
+        let point = BLS12381Curve::generator();
+        let icicle_point = point.to_icicle();
+        let res =
+            ShortWeierstrassProjectivePoint::<BLS12381Curve>::from_icicle_affine(&icicle_point)
+                .unwrap();
+        assert_eq!(point, res)
+    }
+
+    #[test]
+    fn icicle_g1_msm() {
+        const LEN: usize = 20;
+        let eight: BLS12381FieldElement = FieldElement::from(8);
+        let lambda_scalars = vec![eight; LEN];
+        let lambda_points = (0..LEN).map(|_| point_times_5()).collect::<Vec<_>>();
+        let expected = msm(&lambda_scalars, &lambda_points).unwrap();
+        let res = bls12_381_g1_msm(&lambda_scalars, &lambda_points, None).unwrap();
+        assert_eq!(res, expected);
+    }
+}
diff --git a/math/src/gpu/mod.rs b/math/src/gpu/mod.rs
index ee1867d80..a19c098b2 100644
--- a/math/src/gpu/mod.rs
+++ b/math/src/gpu/mod.rs
@@ -6,3 +6,6 @@ If you were using the `--all-features` flag please read this crate's Cargo.toml"
 
 #[cfg(feature = "cuda")]
 pub mod cuda;
+
+#[cfg(feature = "icicle")]
+pub mod icicle;
diff --git a/math/src/msm/naive.rs b/math/src/msm/naive.rs
index fec773073..50a1bad9b 100644
--- a/math/src/msm/naive.rs
+++ b/math/src/msm/naive.rs
@@ -2,16 +2,22 @@ use core::fmt::Display;
 
 use crate::cyclic_group::IsGroup;
 use crate::unsigned_integer::traits::IsUnsignedInteger;
+#[cfg(feature = "icicle")]
+use icicle_core::error::IcicleError;
 
 #[derive(Debug)]
 pub enum MSMError {
     LengthMismatch(usize, usize),
+    #[cfg(feature = "icicle")]
+    Icicle(IcicleError),
 }
 
 impl Display for MSMError {
     fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
         match self {
             MSMError::LengthMismatch(cs, points) => write!(f, "`cs` and `points` must be of the same length to compute `msm`. Got: {cs} and {points}"),
+            #[cfg(feature = "icicle")]
+            MSMError::Icicle(e) => write!(f, "Icicle GPU backend failure. {:?}", e),
         }
     }
 }
diff --git a/math/src/msm/pippenger.rs b/math/src/msm/pippenger.rs
index 1ccfce9c3..8e7c1ba59 100644
--- a/math/src/msm/pippenger.rs
+++ b/math/src/msm/pippenger.rs
@@ -1,7 +1,14 @@
-use crate::{cyclic_group::IsGroup, unsigned_integer::element::UnsignedInteger};
+use crate::{
+    cyclic_group::IsGroup,
+    field::{element::FieldElement, traits::IsField},
+    unsigned_integer::element::UnsignedInteger,
+    gpu::icicle::GpuMSMPoint
+};
 
 use super::naive::MSMError;
-
+#[cfg(feature = "icicle")]
+use crate::gpu::icicle::icicle_msm;
+use crate::traits::ByteConversion;
 use alloc::vec;
 
 /// This function computes the multiscalar multiplication (MSM).
@@ -15,20 +22,39 @@ use alloc::vec;
 /// If `points` and `cs` are empty, then `msm` returns the zero element of the group.
 ///
 /// Panics if `cs` and `points` have different lengths.
-pub fn msm<const NUM_LIMBS: usize, G>(
-    cs: &[UnsignedInteger<NUM_LIMBS>],
+pub fn msm<const NUM_LIMBS: usize, F: IsField<BaseType = UnsignedInteger<NUM_LIMBS>>, G>(
+    cs: &[FieldElement<F>],
     points: &[G],
 ) -> Result<G, MSMError>
 where
-    G: IsGroup,
+    G: IsGroup + GpuMSMPoint,
+    FieldElement<F>: ByteConversion,
 {
     if cs.len() != points.len() {
         return Err(MSMError::LengthMismatch(cs.len(), points.len()));
     }
 
-    let window_size = optimum_window_size(cs.len());
+    #[cfg(feature = "icicle")]
+    {
+        if !G::curve_name().is_empty() {
+            icicle_msm(cs, points)
+        } else {
+            println!(
+                "Icicle msm failed for field {}. Program will fallback to CPU.",
+                core::any::type_name::<F>()
+            );
+            let window_size = optimum_window_size(cs.len());
+            let cs = cs.iter().map(|cs| *cs.value()).collect::<Vec<_>>();
+            Ok(msm_with(&cs, points, window_size))
+        }
+    }
 
-    Ok(msm_with(cs, points, window_size))
+    #[cfg(not(feature = "icicle"))]
+    {
+        let window_size = optimum_window_size(cs.len());
+        let cs = cs.iter().map(|cs| *cs.value()).collect::<Vec<_>>();
+        Ok(msm_with(&cs, points, window_size))
+    }
 }
 
 fn optimum_window_size(data_length: usize) -> usize {
diff --git a/provers/groth16/src/prover.rs b/provers/groth16/src/prover.rs
index b0dbeb579..559bfb5f5 100644
--- a/provers/groth16/src/prover.rs
+++ b/provers/groth16/src/prover.rs
@@ -66,29 +66,20 @@ impl Proof {
 pub struct Prover;
 impl Prover {
     pub fn prove(w: &[FrElement], qap: &QuadraticArithmeticProgram, pk: &ProvingKey) -> Proof {
-        let h_coefficients = qap
-            .calculate_h_coefficients(w)
-            .iter()
-            .map(|elem| elem.representative())
-            .collect::<Vec<_>>();
-
-        let w = w
-            .iter()
-            .map(|elem| elem.representative())
-            .collect::<Vec<_>>();
+        let h_coefficients = qap.calculate_h_coefficients(w);
 
         // Sample randomness for hiding
         let r = sample_fr_elem();
         let s = sample_fr_elem();
 
         // [π_1]_1
-        let pi1 = msm(&w, &pk.l_tau_g1)
+        let pi1 = msm(w, &pk.l_tau_g1)
             .unwrap()
             .operate_with(&pk.alpha_g1)
             .operate_with(&pk.delta_g1.operate_with_self(r.representative()));
 
         // [π_2]_2
-        let pi2 = msm(&w, &pk.r_tau_g2)
+        let pi2 = msm(w, &pk.r_tau_g2)
             .unwrap()
             .operate_with(&pk.beta_g2)
             .operate_with(&pk.delta_g2.operate_with_self(s.representative()));
@@ -108,7 +99,7 @@ impl Prover {
         .unwrap();
 
         // [π_2]_1
-        let pi2_g1 = msm(&w, &pk.r_tau_g1)
+        let pi2_g1 = msm(w, &pk.r_tau_g1)
             .unwrap()
             .operate_with(&pk.beta_g1)
             .operate_with(&pk.delta_g1.operate_with_self(s.representative()));
diff --git a/provers/groth16/src/verifier.rs b/provers/groth16/src/verifier.rs
index b83c4b215..0c0728475 100644
--- a/provers/groth16/src/verifier.rs
+++ b/provers/groth16/src/verifier.rs
@@ -6,14 +6,7 @@ use crate::setup::VerifyingKey;
 
 pub fn verify(vk: &VerifyingKey, proof: &Proof, pub_inputs: &[FrElement]) -> bool {
     // [γ^{-1} * (β*l(τ) + α*r(τ) + o(τ))]_1
-    let k_tau_assigned_verifier_g1 = msm(
-        &pub_inputs
-            .iter()
-            .map(|elem| elem.representative())
-            .collect::<Vec<_>>(),
-        &vk.verifier_k_tau_g1,
-    )
-    .unwrap();
+    let k_tau_assigned_verifier_g1 = msm(pub_inputs, &vk.verifier_k_tau_g1).unwrap();
 
     Pairing::compute(&proof.pi3, &vk.delta_g2).unwrap()
         * vk.alpha_g1_times_beta_g2.clone()
diff --git a/provers/plonk/Cargo.toml b/provers/plonk/Cargo.toml
index 5f3784d48..e66804397 100644
--- a/provers/plonk/Cargo.toml
+++ b/provers/plonk/Cargo.toml
@@ -9,6 +9,7 @@ edition = "2021"
 lambdaworks-math.workspace = true
 lambdaworks-crypto.workspace = true
 serde = { version = "1.0", features = ["derive"] }
+icicle-core = { git = "https://github.com/ingonyama-zk/icicle.git", tag = "v1.4.0" }
 serde_json = "1.0"
 sha3 = { version = "0.10", default-features = false }
 sha2 = { version = "0.10", default-features = false }
diff --git a/provers/plonk/src/prover.rs b/provers/plonk/src/prover.rs
index 3281c48db..aa2654e0c 100644
--- a/provers/plonk/src/prover.rs
+++ b/provers/plonk/src/prover.rs
@@ -15,6 +15,9 @@ use lambdaworks_math::{
 };
 use lambdaworks_math::{field::traits::IsField, traits::ByteConversion};
 
+//TODO: feature gate
+use lambdaworks_math::gpu::icicle::IcicleFFT;
+
 /// Plonk proof.
 /// The challenges are denoted
 ///     Round 2: β,γ,
@@ -279,7 +282,7 @@ struct Round5Result<F: IsField, Hiding> {
 
 impl<F, CS, R> Prover<F, CS, R>
 where
-    F: IsField + IsFFTField,
+    F: IsField + IsFFTField + IcicleFFT,
     CS: IsCommitmentScheme<F>,
     FieldElement<F>: ByteConversion,
     CS::Commitment: AsBytes,
diff --git a/provers/plonk/src/setup.rs b/provers/plonk/src/setup.rs
index 01f08588f..c8398415d 100644
--- a/provers/plonk/src/setup.rs
+++ b/provers/plonk/src/setup.rs
@@ -11,6 +11,9 @@ use lambdaworks_math::field::{element::FieldElement, traits::IsField};
 use lambdaworks_math::polynomial::Polynomial;
 use lambdaworks_math::traits::{AsBytes, ByteConversion};
 
+//TODO: feature gate
+use lambdaworks_math::gpu::icicle::IcicleFFT;
+
 // TODO: implement getters
 pub struct Witness<F: IsField> {
     pub a: Vec<FieldElement<F>>,
@@ -34,7 +37,10 @@ impl<F: IsField> Witness<F> {
 
 // TODO: implement getters
 #[derive(Clone)]
-pub struct CommonPreprocessedInput<F: IsField> {
+pub struct CommonPreprocessedInput<F: IsField + IcicleFFT>
+where
+    FieldElement<F>: ByteConversion,
+{
     pub n: usize,
     /// Number of constraints
     pub domain: Vec<FieldElement<F>>,
@@ -56,7 +62,10 @@ pub struct CommonPreprocessedInput<F: IsField> {
     pub s3_lagrange: Vec<FieldElement<F>>,
 }
 
-impl<F: IsFFTField> CommonPreprocessedInput<F> {
+impl<F: IsFFTField + IcicleFFT> CommonPreprocessedInput<F>
+where
+    FieldElement<F>: ByteConversion,
+{
     pub fn from_constraint_system(
         system: &ConstraintSystem<F>,
         order_r_minus_1_root_unity: &FieldElement<F>,
@@ -113,10 +122,13 @@ pub struct VerificationKey<G1Point> {
     pub s3_1: G1Point,
 }
 
-pub fn setup<F: IsField, CS: IsCommitmentScheme<F>>(
+pub fn setup<F: IsField + IcicleFFT, CS: IsCommitmentScheme<F>>(
     common_input: &CommonPreprocessedInput<F>,
     commitment_scheme: &CS,
-) -> VerificationKey<CS::Commitment> {
+) -> VerificationKey<CS::Commitment>
+where
+    FieldElement<F>: ByteConversion,
+{
     VerificationKey {
         qm_1: commitment_scheme.commit(&common_input.qm),
         ql_1: commitment_scheme.commit(&common_input.ql),
diff --git a/provers/plonk/src/verifier.rs b/provers/plonk/src/verifier.rs
index 8393148dd..cf3bff2e6 100644
--- a/provers/plonk/src/verifier.rs
+++ b/provers/plonk/src/verifier.rs
@@ -9,12 +9,18 @@ use std::marker::PhantomData;
 use crate::prover::Proof;
 use crate::setup::{new_strong_fiat_shamir_transcript, CommonPreprocessedInput, VerificationKey};
 
+//TODO: feature gate
+use lambdaworks_math::gpu::icicle::IcicleFFT;
+
 pub struct Verifier<F: IsField, CS: IsCommitmentScheme<F>> {
     commitment_scheme: CS,
     phantom: PhantomData<F>,
 }
 
-impl<F: IsField + IsFFTField, CS: IsCommitmentScheme<F>> Verifier<F, CS> {
+impl<F: IsField + IsFFTField + IcicleFFT, CS: IsCommitmentScheme<F>> Verifier<F, CS>
+where
+    FieldElement<F>: ByteConversion,
+{
     pub fn new(commitment_scheme: CS) -> Self {
         Self {
             commitment_scheme,
diff --git a/provers/stark/Cargo.toml b/provers/stark/Cargo.toml
index 2c181a729..a8245319c 100644
--- a/provers/stark/Cargo.toml
+++ b/provers/stark/Cargo.toml
@@ -11,6 +11,7 @@ crate-type = ["cdylib", "rlib"]
 [dependencies]
 lambdaworks-math = { workspace = true , features = ["std", "lambdaworks-serde-binary"] }
 lambdaworks-crypto = { workspace = true, features = ["std", "serde"] }
+icicle-core = { git = "https://github.com/ingonyama-zk/icicle.git", tag = "v1.4.0" }
 miden-core = { git="https://github.com/lambdaclass/miden-vm", optional=true}
 
 rand = "0.8.5"
diff --git a/provers/stark/src/constraints/evaluator.rs b/provers/stark/src/constraints/evaluator.rs
index 140a104ab..4a8032b79 100644
--- a/provers/stark/src/constraints/evaluator.rs
+++ b/provers/stark/src/constraints/evaluator.rs
@@ -8,7 +8,12 @@ use crate::{frame::Frame, prover::evaluate_polynomial_on_lde_domain};
 use itertools::Itertools;
 #[cfg(all(debug_assertions, not(feature = "parallel")))]
 use lambdaworks_math::polynomial::Polynomial;
-use lambdaworks_math::{fft::errors::FFTError, field::element::FieldElement, traits::AsBytes};
+use lambdaworks_math::{
+    fft::errors::FFTError,
+    field::element::FieldElement,
+    gpu::icicle::IcicleFFT,
+    traits::{AsBytes, ByteConversion},
+};
 #[cfg(feature = "parallel")]
 use rayon::{
     iter::IndexedParallelIterator,
@@ -39,9 +44,10 @@ impl<A: AIR> ConstraintEvaluator<A> {
         rap_challenges: &[FieldElement<A::FieldExtension>],
     ) -> Vec<FieldElement<A::FieldExtension>>
     where
-        FieldElement<A::Field>: AsBytes + Send + Sync,
+        FieldElement<A::Field>: AsBytes + Send + Sync + ByteConversion,
         FieldElement<A::FieldExtension>: AsBytes + Send + Sync,
         A: Send + Sync,
+        A::Field: IcicleFFT,
     {
         let boundary_constraints = &self.boundary_constraints;
         let number_of_b_constraints = boundary_constraints.constraints.len();
diff --git a/provers/stark/src/constraints/transition.rs b/provers/stark/src/constraints/transition.rs
index 52fe0ba70..09ed5a63a 100644
--- a/provers/stark/src/constraints/transition.rs
+++ b/provers/stark/src/constraints/transition.rs
@@ -4,7 +4,9 @@ use crate::prover::evaluate_polynomial_on_lde_domain;
 use itertools::Itertools;
 use lambdaworks_math::field::element::FieldElement;
 use lambdaworks_math::field::traits::{IsFFTField, IsField, IsSubFieldOf};
+use lambdaworks_math::gpu::icicle::IcicleFFT;
 use lambdaworks_math::polynomial::Polynomial;
+use lambdaworks_math::traits::ByteConversion;
 use num_integer::Integer;
 use std::ops::Div;
 /// TransitionConstraint represents the behaviour that a transition constraint
@@ -105,7 +107,11 @@ where
     }
 
     /// Compute evaluations of the constraints zerofier over a LDE domain.
-    fn zerofier_evaluations_on_extended_domain(&self, domain: &Domain<F>) -> Vec<FieldElement<F>> {
+    fn zerofier_evaluations_on_extended_domain(&self, domain: &Domain<F>) -> Vec<FieldElement<F>>
+    where
+        F: IcicleFFT,
+        FieldElement<F>: ByteConversion,
+    {
         let blowup_factor = domain.blowup_factor;
         let trace_length = domain.trace_roots_of_unity.len();
         let trace_primitive_root = &domain.trace_primitive_root;
diff --git a/provers/stark/src/debug.rs b/provers/stark/src/debug.rs
index fc62257a4..c4e2aa80e 100644
--- a/provers/stark/src/debug.rs
+++ b/provers/stark/src/debug.rs
@@ -6,7 +6,9 @@ use lambdaworks_math::{
         element::FieldElement,
         traits::{IsFFTField, IsField},
     },
+    gpu::icicle::IcicleFFT,
     polynomial::Polynomial,
+    traits::ByteConversion,
 };
 use log::{error, info};
 
@@ -17,7 +19,13 @@ pub fn validate_trace<A: AIR>(
     aux_trace_polys: &[Polynomial<FieldElement<A::FieldExtension>>],
     domain: &Domain<A::Field>,
     rap_challenges: &[FieldElement<A::FieldExtension>],
-) -> bool {
+) -> bool
+where
+    <A as AIR>::Field: IcicleFFT,
+    FieldElement<<A as AIR>::Field>: ByteConversion,
+    FieldElement<<A as AIR>::FieldExtension>: ByteConversion,
+    <A as AIR>::FieldExtension: IsFFTField + IcicleFFT,
+{
     info!("Starting constraints validation over trace...");
     let mut ret = true;
 
diff --git a/provers/stark/src/examples/dummy_air.rs b/provers/stark/src/examples/dummy_air.rs
index 1206ed153..8102addb9 100644
--- a/provers/stark/src/examples/dummy_air.rs
+++ b/provers/stark/src/examples/dummy_air.rs
@@ -15,6 +15,8 @@ use lambdaworks_math::field::{
     element::FieldElement, fields::fft_friendly::stark_252_prime_field::Stark252PrimeField,
     traits::IsFFTField,
 };
+use lambdaworks_math::gpu::icicle::IcicleFFT;
+use lambdaworks_math::traits::ByteConversion;
 
 type StarkField = Stark252PrimeField;
 
@@ -32,7 +34,8 @@ impl<F: IsFFTField> FibConstraint<F> {
 
 impl<F> TransitionConstraint<F, F> for FibConstraint<F>
 where
-    F: IsFFTField + Send + Sync,
+    F: IsFFTField + Send + Sync + IcicleFFT,
+    FieldElement<F>: ByteConversion,
 {
     fn degree(&self) -> usize {
         1
@@ -81,7 +84,8 @@ impl<F: IsFFTField> BitConstraint<F> {
 
 impl<F> TransitionConstraint<F, F> for BitConstraint<F>
 where
-    F: IsFFTField + Send + Sync,
+    F: IsFFTField + Send + Sync + IcicleFFT,
+    FieldElement<F>: ByteConversion,
 {
     fn degree(&self) -> usize {
         2
diff --git a/provers/stark/src/examples/fibonacci_2_cols_shifted.rs b/provers/stark/src/examples/fibonacci_2_cols_shifted.rs
index 6d3eb76d1..6315b36dd 100644
--- a/provers/stark/src/examples/fibonacci_2_cols_shifted.rs
+++ b/provers/stark/src/examples/fibonacci_2_cols_shifted.rs
@@ -11,7 +11,8 @@ use crate::{
 };
 use lambdaworks_math::{
     field::{element::FieldElement, traits::IsFFTField},
-    traits::AsBytes,
+    gpu::icicle::IcicleFFT,
+    traits::{AsBytes, ByteConversion},
 };
 use std::marker::PhantomData;
 
@@ -30,7 +31,8 @@ impl<F: IsFFTField> ShiftedFibTransition1<F> {
 
 impl<F> TransitionConstraint<F, F> for ShiftedFibTransition1<F>
 where
-    F: IsFFTField + Send + Sync,
+    F: IsFFTField + Send + Sync + IcicleFFT,
+    FieldElement<F>: ByteConversion,
 {
     fn degree(&self) -> usize {
         1
@@ -78,7 +80,8 @@ impl<F: IsFFTField> ShiftedFibTransition2<F> {
 
 impl<F> TransitionConstraint<F, F> for ShiftedFibTransition2<F>
 where
-    F: IsFFTField + Send + Sync,
+    F: IsFFTField + Send + Sync + IcicleFFT,
+    FieldElement<F>: ByteConversion,
 {
     fn degree(&self) -> usize {
         1
@@ -149,7 +152,8 @@ where
 /// for all `i`. Also, `Col0_0` is constrained to be `1`.
 impl<F> AIR for Fibonacci2ColsShifted<F>
 where
-    F: IsFFTField + Send + Sync + 'static,
+    F: IsFFTField + Send + Sync + 'static + IcicleFFT,
+    FieldElement<F>: ByteConversion,
 {
     type Field = F;
     type FieldExtension = F;
diff --git a/provers/stark/src/examples/fibonacci_2_columns.rs b/provers/stark/src/examples/fibonacci_2_columns.rs
index 0cd2a789a..9846cecb4 100644
--- a/provers/stark/src/examples/fibonacci_2_columns.rs
+++ b/provers/stark/src/examples/fibonacci_2_columns.rs
@@ -12,7 +12,11 @@ use crate::{
     trace::TraceTable,
     traits::AIR,
 };
-use lambdaworks_math::field::{element::FieldElement, traits::IsFFTField};
+use lambdaworks_math::{
+    field::{element::FieldElement, traits::IsFFTField},
+    gpu::icicle::IcicleFFT,
+    traits::ByteConversion,
+};
 
 #[derive(Clone)]
 struct FibTransition1<F: IsFFTField> {
@@ -29,7 +33,8 @@ impl<F: IsFFTField> FibTransition1<F> {
 
 impl<F> TransitionConstraint<F, F> for FibTransition1<F>
 where
-    F: IsFFTField + Send + Sync,
+    F: IsFFTField + Send + Sync + IcicleFFT,
+    FieldElement<F>: ByteConversion,
 {
     fn degree(&self) -> usize {
         1
@@ -79,7 +84,8 @@ impl<F: IsFFTField> FibTransition2<F> {
 
 impl<F> TransitionConstraint<F, F> for FibTransition2<F>
 where
-    F: IsFFTField + Send + Sync,
+    F: IsFFTField + Send + Sync + IcicleFFT,
+    FieldElement<F>: ByteConversion,
 {
     fn degree(&self) -> usize {
         1
@@ -128,7 +134,8 @@ where
 /// stacked in row-major order.
 impl<F> AIR for Fibonacci2ColsAIR<F>
 where
-    F: IsFFTField + Send + Sync + 'static,
+    F: IsFFTField + Send + Sync + 'static + IcicleFFT,
+    FieldElement<F>: ByteConversion,
 {
     type Field = F;
     type FieldExtension = F;
diff --git a/provers/stark/src/examples/fibonacci_rap.rs b/provers/stark/src/examples/fibonacci_rap.rs
index 18084f86f..d2a1cd7d3 100644
--- a/provers/stark/src/examples/fibonacci_rap.rs
+++ b/provers/stark/src/examples/fibonacci_rap.rs
@@ -14,6 +14,7 @@ use crate::{
 use lambdaworks_crypto::fiat_shamir::is_transcript::IsTranscript;
 use lambdaworks_math::{
     field::{element::FieldElement, traits::IsFFTField},
+    gpu::icicle::IcicleFFT,
     helpers::resize_to_next_power_of_two,
     traits::ByteConversion,
 };
@@ -33,7 +34,8 @@ impl<F: IsFFTField> FibConstraint<F> {
 
 impl<F> TransitionConstraint<F, F> for FibConstraint<F>
 where
-    F: IsFFTField + Send + Sync,
+    F: IsFFTField + Send + Sync + IcicleFFT,
+    FieldElement<F>: ByteConversion,
 {
     fn degree(&self) -> usize {
         1
@@ -85,7 +87,8 @@ impl<F: IsFFTField> PermutationConstraint<F> {
 
 impl<F> TransitionConstraint<F, F> for PermutationConstraint<F>
 where
-    F: IsFFTField + Send + Sync,
+    F: IsFFTField + Send + Sync + IcicleFFT,
+    FieldElement<F>: ByteConversion,
 {
     fn degree(&self) -> usize {
         2
@@ -145,7 +148,7 @@ where
 
 impl<F> AIR for FibonacciRAP<F>
 where
-    F: IsFFTField + Send + Sync + 'static,
+    F: IsFFTField + Send + Sync + 'static + IcicleFFT,
     FieldElement<F>: ByteConversion,
 {
     type Field = F;
diff --git a/provers/stark/src/examples/quadratic_air.rs b/provers/stark/src/examples/quadratic_air.rs
index 3d07d78e3..f64b35400 100644
--- a/provers/stark/src/examples/quadratic_air.rs
+++ b/provers/stark/src/examples/quadratic_air.rs
@@ -11,7 +11,11 @@ use crate::{
     trace::TraceTable,
     traits::AIR,
 };
-use lambdaworks_math::field::{element::FieldElement, traits::IsFFTField};
+use lambdaworks_math::{
+    field::{element::FieldElement, traits::IsFFTField},
+    gpu::icicle::IcicleFFT,
+    traits::ByteConversion,
+};
 
 #[derive(Clone)]
 struct QuadraticConstraint<F: IsFFTField> {
@@ -28,7 +32,8 @@ impl<F: IsFFTField> QuadraticConstraint<F> {
 
 impl<F> TransitionConstraint<F, F> for QuadraticConstraint<F>
 where
-    F: IsFFTField + Send + Sync,
+    F: IsFFTField + Send + Sync + IcicleFFT,
+    FieldElement<F>: ByteConversion,
 {
     fn degree(&self) -> usize {
         2
@@ -81,7 +86,8 @@ where
 
 impl<F> AIR for QuadraticAIR<F>
 where
-    F: IsFFTField + Send + Sync + 'static,
+    F: IsFFTField + Send + Sync + 'static + IcicleFFT,
+    FieldElement<F>: ByteConversion,
 {
     type Field = F;
     type FieldExtension = F;
diff --git a/provers/stark/src/examples/simple_fibonacci.rs b/provers/stark/src/examples/simple_fibonacci.rs
index 204aa938c..72b490ea0 100644
--- a/provers/stark/src/examples/simple_fibonacci.rs
+++ b/provers/stark/src/examples/simple_fibonacci.rs
@@ -9,7 +9,11 @@ use crate::{
     trace::TraceTable,
     traits::AIR,
 };
-use lambdaworks_math::field::{element::FieldElement, traits::IsFFTField};
+use lambdaworks_math::{
+    field::{element::FieldElement, traits::IsFFTField},
+    gpu::icicle::IcicleFFT,
+    traits::ByteConversion,
+};
 use std::marker::PhantomData;
 
 #[derive(Clone)]
@@ -27,7 +31,8 @@ impl<F: IsFFTField> FibConstraint<F> {
 
 impl<F> TransitionConstraint<F, F> for FibConstraint<F>
 where
-    F: IsFFTField + Send + Sync,
+    F: IsFFTField + Send + Sync + IcicleFFT,
+    FieldElement<F>: ByteConversion,
 {
     fn degree(&self) -> usize {
         1
@@ -83,7 +88,8 @@ where
 
 impl<F> AIR for FibonacciAIR<F>
 where
-    F: IsFFTField + Send + Sync + 'static,
+    F: IsFFTField + Send + Sync + 'static + IcicleFFT,
+    FieldElement<F>: ByteConversion,
 {
     type Field = F;
     type FieldExtension = F;
diff --git a/provers/stark/src/examples/simple_periodic_cols.rs b/provers/stark/src/examples/simple_periodic_cols.rs
index 369ff2dd0..23c8b05b5 100644
--- a/provers/stark/src/examples/simple_periodic_cols.rs
+++ b/provers/stark/src/examples/simple_periodic_cols.rs
@@ -11,7 +11,11 @@ use crate::{
     trace::TraceTable,
     traits::AIR,
 };
-use lambdaworks_math::field::{element::FieldElement, traits::IsFFTField};
+use lambdaworks_math::{
+    field::{element::FieldElement, traits::IsFFTField},
+    gpu::icicle::IcicleFFT,
+    traits::ByteConversion,
+};
 
 pub struct PeriodicConstraint<F: IsFFTField> {
     phantom: PhantomData<F>,
@@ -31,7 +35,8 @@ impl<F: IsFFTField> Default for PeriodicConstraint<F> {
 
 impl<F> TransitionConstraint<F, F> for PeriodicConstraint<F>
 where
-    F: IsFFTField + Send + Sync,
+    F: IsFFTField + Send + Sync + IcicleFFT,
+    FieldElement<F>: ByteConversion,
 {
     fn degree(&self) -> usize {
         1
@@ -101,7 +106,8 @@ where
 
 impl<F> AIR for SimplePeriodicAIR<F>
 where
-    F: IsFFTField + Send + Sync + 'static,
+    F: IsFFTField + Send + Sync + 'static + IcicleFFT,
+    FieldElement<F>: ByteConversion,
 {
     type Field = F;
     type FieldExtension = F;
diff --git a/provers/stark/src/fri/mod.rs b/provers/stark/src/fri/mod.rs
index a841a8aea..4598b92aa 100644
--- a/provers/stark/src/fri/mod.rs
+++ b/provers/stark/src/fri/mod.rs
@@ -10,7 +10,9 @@ use lambdaworks_math::{
 };
 pub use lambdaworks_math::{
     field::{element::FieldElement, fields::u64_prime_field::U64PrimeField},
+    gpu::icicle::IcicleFFT,
     polynomial::Polynomial,
+    traits::ByteConversion,
 };
 
 use crate::config::{BatchedMerkleTree, BatchedMerkleTreeBackend};
@@ -19,7 +21,7 @@ use self::fri_commitment::FriLayer;
 use self::fri_decommit::FriDecommitment;
 use self::fri_functions::fold_polynomial;
 
-pub fn commit_phase<F: IsFFTField + IsSubFieldOf<E>, E: IsField>(
+pub fn commit_phase<F: IsFFTField + IsSubFieldOf<E>, E: IsField + IsFFTField + IcicleFFT>(
     number_layers: usize,
     p_0: Polynomial<FieldElement<E>>,
     transcript: &mut impl IsTranscript<E>,
@@ -31,7 +33,7 @@ pub fn commit_phase<F: IsFFTField + IsSubFieldOf<E>, E: IsField>(
 )
 where
     FieldElement<F>: AsBytes + Sync + Send,
-    FieldElement<E>: AsBytes + Sync + Send,
+    FieldElement<E>: AsBytes + Sync + Send + ByteConversion,
 {
     let mut domain_size = domain_size;
 
@@ -112,14 +114,14 @@ where
     }
 }
 
-pub fn new_fri_layer<F: IsFFTField + IsSubFieldOf<E>, E: IsField>(
+pub fn new_fri_layer<F: IsFFTField + IsSubFieldOf<E>, E: IsField + IsFFTField + IcicleFFT>(
     poly: &Polynomial<FieldElement<E>>,
     coset_offset: &FieldElement<F>,
     domain_size: usize,
 ) -> crate::fri::fri_commitment::FriLayer<E, BatchedMerkleTreeBackend<E>>
 where
     FieldElement<F>: AsBytes + Sync + Send,
-    FieldElement<E>: AsBytes + Sync + Send,
+    FieldElement<E>: AsBytes + Sync + Send + ByteConversion,
 {
     let mut evaluation =
         Polynomial::evaluate_offset_fft(poly, 1, Some(domain_size), coset_offset).unwrap(); // TODO: return error
diff --git a/provers/stark/src/prover.rs b/provers/stark/src/prover.rs
index 7a5a4210e..746e9a41c 100644
--- a/provers/stark/src/prover.rs
+++ b/provers/stark/src/prover.rs
@@ -10,7 +10,9 @@ use lambdaworks_math::field::traits::{IsField, IsSubFieldOf};
 use lambdaworks_math::traits::AsBytes;
 use lambdaworks_math::{
     field::{element::FieldElement, traits::IsFFTField},
+    gpu::icicle::IcicleFFT,
     polynomial::Polynomial,
+    traits::ByteConversion,
 };
 use log::info;
 
@@ -154,7 +156,8 @@ pub fn evaluate_polynomial_on_lde_domain<F, E>(
 ) -> Result<Vec<FieldElement<E>>, FFTError>
 where
     F: IsFFTField + IsSubFieldOf<E>,
-    E: IsField,
+    E: IsField + IsFFTField + IcicleFFT,
+    FieldElement<E>: ByteConversion,
 {
     let evaluations = Polynomial::evaluate_offset_fft(p, blowup_factor, Some(domain_size), offset)?;
     let step = evaluations.len() / (domain_size * blowup_factor);
@@ -203,9 +206,9 @@ pub trait IsStarkProver<A: AIR> {
     )
     where
         FieldElement<A::Field>: AsBytes + Send + Sync,
-        FieldElement<E>: AsBytes + Send + Sync,
+        FieldElement<E>: AsBytes + Send + Sync + ByteConversion,
         FieldElement<A::FieldExtension>: AsBytes + Send + Sync,
-        E: IsSubFieldOf<A::FieldExtension>,
+        E: IsSubFieldOf<A::FieldExtension> + IcicleFFT + IsFFTField,
         A::Field: IsSubFieldOf<E>,
     {
         // Interpolate columns of `trace`.
@@ -243,8 +246,8 @@ pub trait IsStarkProver<A: AIR> {
     ) -> Vec<Vec<FieldElement<E>>>
     where
         FieldElement<A::Field>: Send + Sync,
-        FieldElement<E>: Send + Sync,
-        E: IsSubFieldOf<A::FieldExtension>,
+        FieldElement<E>: Send + Sync + ByteConversion,
+        E: IsSubFieldOf<A::FieldExtension> + IsFFTField + IcicleFFT,
         A::Field: IsSubFieldOf<E>,
     {
         #[cfg(not(feature = "parallel"))]
@@ -273,8 +276,10 @@ pub trait IsStarkProver<A: AIR> {
         transcript: &mut impl IsTranscript<A::FieldExtension>,
     ) -> Result<Round1<A>, ProvingError>
     where
-        FieldElement<A::Field>: AsBytes + Send + Sync,
-        FieldElement<A::FieldExtension>: AsBytes + Send + Sync,
+        A::Field: IcicleFFT,
+        A::FieldExtension: IsFFTField + IcicleFFT,
+        FieldElement<A::Field>: AsBytes + Send + Sync + ByteConversion,
+        FieldElement<A::FieldExtension>: AsBytes + Send + Sync + ByteConversion,
     {
         let (trace_polys, evaluations, main_merkle_tree, main_merkle_root) =
             Self::interpolate_and_commit::<A::Field>(main_trace, domain, transcript);
@@ -358,8 +363,10 @@ pub trait IsStarkProver<A: AIR> {
     ) -> Round2<A::FieldExtension>
     where
         A: Send + Sync,
-        FieldElement<A::Field>: AsBytes + Send + Sync,
-        FieldElement<A::FieldExtension>: AsBytes + Send + Sync,
+        A::Field: IcicleFFT,
+        A::FieldExtension: IsFFTField + IcicleFFT,
+        FieldElement<A::Field>: AsBytes + Send + Sync + ByteConversion,
+        FieldElement<A::FieldExtension>: AsBytes + Send + Sync + ByteConversion,
     {
         // Compute the evaluations of the composition polynomial on the LDE domain.
         let evaluator = ConstraintEvaluator::new(air, &round_1_result.rap_challenges);
@@ -464,8 +471,9 @@ pub trait IsStarkProver<A: AIR> {
         transcript: &mut impl IsTranscript<A::FieldExtension>,
     ) -> Round4<A::Field, A::FieldExtension>
     where
+        A::FieldExtension: IsFFTField + IcicleFFT,
         FieldElement<A::Field>: AsBytes + Send + Sync,
-        FieldElement<A::FieldExtension>: AsBytes + Send + Sync,
+        FieldElement<A::FieldExtension>: AsBytes + Send + Sync + ByteConversion,
     {
         let coset_offset_u64 = air.context().proof_options.coset_offset;
         let coset_offset = FieldElement::<A::Field>::from(coset_offset_u64);
@@ -797,8 +805,10 @@ pub trait IsStarkProver<A: AIR> {
     ) -> Result<StarkProof<A::Field, A::FieldExtension>, ProvingError>
     where
         A: Send + Sync,
-        FieldElement<A::Field>: AsBytes + Send + Sync,
-        FieldElement<A::FieldExtension>: AsBytes + Send + Sync,
+        A::Field: IcicleFFT,
+        A::FieldExtension: IsFFTField + IcicleFFT,
+        FieldElement<A::Field>: AsBytes + Send + Sync + ByteConversion,
+        FieldElement<A::FieldExtension>: AsBytes + Send + Sync + ByteConversion,
     {
         info!("Started proof generation...");
         #[cfg(feature = "instruments")]
diff --git a/provers/stark/src/trace.rs b/provers/stark/src/trace.rs
index b9e3de0cf..9093fa78f 100644
--- a/provers/stark/src/trace.rs
+++ b/provers/stark/src/trace.rs
@@ -4,7 +4,9 @@ use lambdaworks_math::fft::errors::FFTError;
 use lambdaworks_math::field::traits::{IsField, IsSubFieldOf};
 use lambdaworks_math::{
     field::{element::FieldElement, traits::IsFFTField},
+    gpu::icicle::IcicleFFT,
     polynomial::Polynomial,
+    traits::ByteConversion,
 };
 #[cfg(feature = "parallel")]
 use rayon::prelude::{IntoParallelRefIterator, ParallelIterator};
@@ -137,7 +139,8 @@ impl<F: IsField> TraceTable<F> {
     pub fn compute_trace_polys<S>(&self) -> Vec<Polynomial<FieldElement<F>>>
     where
         S: IsFFTField + IsSubFieldOf<F>,
-        FieldElement<F>: Send + Sync,
+        FieldElement<F>: Send + Sync + ByteConversion,
+        F: IcicleFFT + IsFFTField,
     {
         let columns = self.columns();
         #[cfg(feature = "parallel")]
diff --git a/provers/stark/src/traits.rs b/provers/stark/src/traits.rs
index d33343bb2..f78c1e425 100644
--- a/provers/stark/src/traits.rs
+++ b/provers/stark/src/traits.rs
@@ -6,7 +6,9 @@ use lambdaworks_math::{
         element::FieldElement,
         traits::{IsFFTField, IsField, IsSubFieldOf},
     },
+    gpu::icicle::IcicleFFT,
     polynomial::Polynomial,
+    traits::ByteConversion,
 };
 
 use crate::{constraints::transition::TransitionConstraint, domain::Domain};
@@ -124,7 +126,11 @@ pub trait AIR {
         vec![]
     }
 
-    fn get_periodic_column_polynomials(&self) -> Vec<Polynomial<FieldElement<Self::Field>>> {
+    fn get_periodic_column_polynomials(&self) -> Vec<Polynomial<FieldElement<Self::Field>>>
+    where
+        <Self as AIR>::Field: IcicleFFT,
+        FieldElement<<Self as AIR>::Field>: ByteConversion,
+    {
         let mut result = Vec::new();
         for periodic_column in self.get_periodic_column_values() {
             let values: Vec<_> = periodic_column
@@ -150,7 +156,11 @@ pub trait AIR {
     fn transition_zerofier_evaluations(
         &self,
         domain: &Domain<Self::Field>,
-    ) -> Vec<Vec<FieldElement<Self::Field>>> {
+    ) -> Vec<Vec<FieldElement<Self::Field>>>
+    where
+        Self::Field: IcicleFFT,
+        FieldElement<Self::Field>: ByteConversion,
+    {
         let mut evals = vec![Vec::new(); self.num_transition_constraints()];
 
         let mut zerofier_groups: HashMap<ZerofierGroupKey, Vec<FieldElement<Self::Field>>> =
diff --git a/provers/stark/src/verifier.rs b/provers/stark/src/verifier.rs
index ea25dbe79..d4afa89d9 100644
--- a/provers/stark/src/verifier.rs
+++ b/provers/stark/src/verifier.rs
@@ -14,7 +14,8 @@ use lambdaworks_math::{
         element::FieldElement,
         traits::{IsFFTField, IsField, IsSubFieldOf},
     },
-    traits::AsBytes,
+    gpu::icicle::IcicleFFT,
+    traits::{AsBytes, ByteConversion},
 };
 #[cfg(not(feature = "test_fiat_shamir"))]
 use log::error;
@@ -216,7 +217,11 @@ pub trait IsStarkVerifier<A: AIR> {
         proof: &StarkProof<A::Field, A::FieldExtension>,
         domain: &Domain<A::Field>,
         challenges: &Challenges<A>,
-    ) -> bool {
+    ) -> bool
+    where
+        A::Field: IcicleFFT,
+        FieldElement<<A as AIR>::Field>: ByteConversion,
+    {
         let boundary_constraints = air.boundary_constraints(&challenges.rap_challenges);
 
         let trace_length = air.trace_length();
@@ -710,8 +715,9 @@ pub trait IsStarkVerifier<A: AIR> {
         mut transcript: impl IsTranscript<A::FieldExtension>,
     ) -> bool
     where
-        FieldElement<A::Field>: AsBytes + Sync + Send,
+        FieldElement<A::Field>: AsBytes + Sync + Send + ByteConversion,
         FieldElement<A::FieldExtension>: AsBytes + Sync + Send,
+        A::Field: IcicleFFT,
     {
         // Verify there are enough queries
         if proof.query_list.len() < proof_options.fri_number_of_queries {