feat(Field): Add mersenne31

fuzz/no_gpu_fuzz/Cargo.toml

@@ -17,12 +17,21 @@ ibig = "0.3.6"
 p3-goldilocks = { git = "https://github.com/Plonky3/Plonky3", rev = "41cd843" }
 p3-field = { git = "https://github.com/Plonky3/Plonky3", rev = "41cd843" }
 
+p3-mersenne-31 = { git = "https://github.com/Plonky3/Plonky3", rev = "41cd843" }
+p3-field = { git = "https://github.com/Plonky3/Plonky3", rev = "41cd843" }
+
 [[bin]]
 name = "field_fuzzer"
 path = "fuzz_targets/field_fuzzer.rs"
 test = false
 doc = false
 
+[[bin]]
+name = "field_fuzz_mersenne31"
+path = "fuzz_targets/field_mersenne31.rs"
+test = false
+doc = false
+
 [[bin]]
 name = "field_mini_goldilocks"
 path = "fuzz_targets/field_mini_goldilocks.rs"

fuzz/no_gpu_fuzz/fuzz_targets/field_mersenne31.rs

@@ -0,0 +1,90 @@
+#![no_main]
+
+use libfuzzer_sys::fuzz_target;
+use lambdaworks_math::field::{
+    element::FieldElement, 
+    fields::{
+        mersenne31::field::{Mersenne31Field, MERSENNE_31_PRIME_FIELD_ORDER},
+    }
+};
+use p3_mersenne_31::Mersenne31;
+use p3_field::{Field, PrimeField32, PrimeField64, AbstractField};
+
+fuzz_target!(|values: (u32, u32)| {
+    // Note: we filter values outside of order as it triggers an assert within plonky3 disallowing values n >= Self::Order
+    if values.0 >= MERSENNE_31_PRIME_FIELD_ORDER || values.1 >= MERSENNE_31_PRIME_FIELD_ORDER {
+        return
+    }
+
+    let (value_u32_a, value_u32_b) = values;
+
+    let a =  FieldElement::<Mersenne31Field>::from(value_u32_a as u64);
+    let b =  FieldElement::<Mersenne31Field>::from(value_u32_b as u64);
+
+    // Note: if we parse using from_canonical_u32 fails due to check that n < Self::Order
+    let a_expected = Mersenne31::from_canonical_u32(value_u32_a);
+    let b_expected = Mersenne31::from_canonical_u32(value_u32_b);
+
+    let add_u32 = &a + &b;
+    let addition = a_expected + b_expected;
+
+    assert_eq!(add_u32.representative(), addition.as_canonical_u32());
+
+    let sub_u32 = &a - &b;
+    let substraction = a_expected - b_expected;
+    assert_eq!(sub_u32.representative(), substraction.as_canonical_u32());
+
+    let mul_u32 = &a * &b;
+    let multiplication = a_expected * b_expected;
+    assert_eq!(mul_u32.representative(), multiplication.as_canonical_u32());
+
+    let pow = &a.pow(b.representative());
+    let expected_pow = a_expected.exp_u64(b_expected.as_canonical_u64());
+    assert_eq!(pow.representative(), expected_pow.as_canonical_u32());
+
+    if value_u32_b != 0 && b.inv().is_ok() && b_expected.try_inverse().is_some() {
+        let div = &a / &b; 
+        assert_eq!(&div * &b, a.clone());
+        let expected_div = a_expected / b_expected;
+        assert_eq!(div.representative(), expected_div.as_canonical_u32());
+    }
+
+    for n in [&a, &b] {
+        match n.sqrt() {
+            Some((fst_sqrt, snd_sqrt)) => {
+                assert_eq!(fst_sqrt.square(), snd_sqrt.square(), "Squared roots don't match each other");
+                assert_eq!(n, &fst_sqrt.square(), "Squared roots don't match original number");
+            }
+            None => {}
+        };
+    }
+
+    // Axioms soundness
+
+    let one = FieldElement::<Mersenne31Field>::one();
+    let zero = FieldElement::<Mersenne31Field>::zero();
+
+    assert_eq!(&a + &zero, a, "Neutral add element a failed");
+    assert_eq!(&b + &zero, b, "Neutral mul element b failed");
+    assert_eq!(&a * &one, a, "Neutral add element a failed");
+    assert_eq!(&b * &one, b, "Neutral mul element b failed");
+
+    assert_eq!(&a + &b, &b + &a, "Commutative add property failed");
+    assert_eq!(&a * &b, &b * &a, "Commutative mul property failed");
+
+    let c = &a * &b;
+    assert_eq!((&a + &b) + &c, &a + (&b + &c), "Associative add property failed");
+    assert_eq!((&a * &b) * &c, &a * (&b * &c), "Associative mul property failed");
+
+    assert_eq!(&a * (&b + &c), &a * &b + &a * &c, "Distributive property failed");
+
+    assert_eq!(&a - &a, zero, "Inverse add a failed");
+    assert_eq!(&b - &b, zero, "Inverse add b failed");
+
+    if a != zero {
+        assert_eq!(&a * a.inv().unwrap(), one, "Inverse mul a failed");
+    }
+    if b != zero {
+        assert_eq!(&b * b.inv().unwrap(), one, "Inverse mul b failed");
+    }
+});

math/benches/criterion_field.rs

@@ -2,6 +2,8 @@ use criterion::{criterion_group, criterion_main, Criterion};
 use pprof::criterion::{Output, PProfProfiler};
 
 mod fields;
+use fields::mersenne31::mersenne31_ops_benchmarks;
+use fields::mersenne31_montgomery::mersenne31_mont_ops_benchmarks;
 use fields::{
     stark252::starkfield_ops_benchmarks, u64_goldilocks::u64_goldilocks_ops_benchmarks,
     u64_goldilocks_montgomery::u64_goldilocks_montgomery_ops_benchmarks,
@@ -10,6 +12,6 @@ use fields::{
 criterion_group!(
     name = field_benches;
     config = Criterion::default().with_profiler(PProfProfiler::new(100, Output::Flamegraph(None)));
-    targets = starkfield_ops_benchmarks, u64_goldilocks_ops_benchmarks, u64_goldilocks_montgomery_ops_benchmarks
+    targets = starkfield_ops_benchmarks, mersenne31_ops_benchmarks, mersenne31_mont_ops_benchmarks, u64_goldilocks_ops_benchmarks, u64_goldilocks_montgomery_ops_benchmarks
 );
 criterion_main!(field_benches);

math/benches/fields/mersenne31.rs

@@ -0,0 +1,195 @@
+use std::hint::black_box;
+
+use criterion::Criterion;
+use lambdaworks_math::field::{element::FieldElement, fields::mersenne31::field::Mersenne31Field};
+use rand::random;
+
+pub type F = FieldElement<Mersenne31Field>;
+
+#[inline(never)]
+#[no_mangle]
+#[export_name = "util::rand_mersenne31_field_elements"]
+pub fn rand_field_elements(num: usize) -> Vec<(F, F)> {
+    let mut result = Vec::with_capacity(num);
+    for _ in 0..result.capacity() {
+        result.push((F::new(random()), F::new(random())));
+    }
+    result
+}
+
+pub fn mersenne31_ops_benchmarks(c: &mut Criterion) {
+    let input: Vec<Vec<(F, F)>> = [1, 10, 100, 1000, 10000, 100000, 1000000]
+        .into_iter()
+        .map(rand_field_elements)
+        .collect::<Vec<_>>();
+    let mut group = c.benchmark_group("Mersenne31 operations");
+
+    for i in input.clone().into_iter() {
+        group.bench_with_input(format!("add {:?}", &i.len()), &i, |bench, i| {
+            bench.iter(|| {
+                for (x, y) in i {
+                    black_box(black_box(x) + black_box(y));
+                }
+            });
+        });
+    }
+
+    for i in input.clone().into_iter() {
+        group.bench_with_input(format!("mul {:?}", &i.len()), &i, |bench, i| {
+            bench.iter(|| {
+                for (x, y) in i {
+                    black_box(black_box(x) * black_box(y));
+                }
+            });
+        });
+    }
+
+    for i in input.clone().into_iter() {
+        group.bench_with_input(format!("pow by 1 {:?}", &i.len()), &i, |bench, i| {
+            bench.iter(|| {
+                for (x, _) in i {
+                    black_box(black_box(x).pow(1_u64));
+                }
+            });
+        });
+    }
+
+    for i in input.clone().into_iter() {
+        group.bench_with_input(format!("square {:?}", &i.len()), &i, |bench, i| {
+            bench.iter(|| {
+                for (x, _) in i {
+                    black_box(black_box(x).square());
+                }
+            });
+        });
+    }
+
+    for i in input.clone().into_iter() {
+        group.bench_with_input(format!("square with pow {:?}", &i.len()), &i, |bench, i| {
+            bench.iter(|| {
+                for (x, _) in i {
+                    black_box(black_box(x).pow(2_u64));
+                }
+            });
+        });
+    }
+
+    for i in input.clone().into_iter() {
+        group.bench_with_input(format!("square with mul {:?}", &i.len()), &i, |bench, i| {
+            bench.iter(|| {
+                for (x, _) in i {
+                    black_box(black_box(x) * black_box(x));
+                }
+            });
+        });
+    }
+
+    for i in input.clone().into_iter() {
+        group.bench_with_input(
+            format!("pow {:?}", &i.len()),
+            &(i, 5u64),
+            |bench, (i, a)| {
+                bench.iter(|| {
+                    for (x, _) in i {
+                        black_box(black_box(x).pow(*a));
+                    }
+                });
+            },
+        );
+    }
+
+    for i in input.clone().into_iter() {
+        group.bench_with_input(format!("sub {:?}", &i.len()), &i, |bench, i| {
+            bench.iter(|| {
+                for (x, y) in i {
+                    black_box(black_box(x) - black_box(y));
+                }
+            });
+        });
+    }
+
+    for i in input.clone().into_iter() {
+        group.bench_with_input(format!("inv {:?}", &i.len()), &i, |bench, i| {
+            bench.iter(|| {
+                for (x, _) in i {
+                    black_box(black_box(x).inv().unwrap());
+                }
+            });
+        });
+    }
+
+    for i in input.clone().into_iter() {
+        group.bench_with_input(format!("div {:?}", &i.len()), &i, |bench, i| {
+            bench.iter(|| {
+                for (x, y) in i {
+                    black_box(black_box(x) / black_box(y));
+                }
+            });
+        });
+    }
+
+    for i in input.clone().into_iter() {
+        group.bench_with_input(format!("eq {:?}", &i.len()), &i, |bench, i| {
+            bench.iter(|| {
+                for (x, y) in i {
+                    black_box(black_box(x) == black_box(y));
+                }
+            });
+        });
+    }
+
+    for i in input.clone().into_iter() {
+        group.bench_with_input(format!("sqrt {:?}", &i.len()), &i, |bench, i| {
+            bench.iter(|| {
+                for (x, _) in i {
+                    black_box(black_box(x).sqrt());
+                }
+            });
+        });
+    }
+
+    for i in input.clone().into_iter() {
+        group.bench_with_input(format!("sqrt squared {:?}", &i.len()), &i, |bench, i| {
+            let i: Vec<F> = i.iter().map(|(x, _)| x * x).collect();
+            bench.iter(|| {
+                for x in &i {
+                    black_box(black_box(x).sqrt());
+                }
+            });
+        });
+    }
+
+    for i in input.clone().into_iter() {
+        group.bench_with_input(format!("bitand {:?}", &i.len()), &i, |bench, i| {
+            // Note: we should strive to have the number of limbs be generic... ideally this benchmark group itself should have a generic type that we call into from the main runner.
+            let i: Vec<(u32, u32)> = i.iter().map(|(x, y)| (*x.value(), *y.value())).collect();
+            bench.iter(|| {
+                for (x, y) in &i {
+                    black_box(black_box(*x) & black_box(*y));
+                }
+            });
+        });
+    }
+
+    for i in input.clone().into_iter() {
+        group.bench_with_input(format!("bitor {:?}", &i.len()), &i, |bench, i| {
+            let i: Vec<(u32, u32)> = i.iter().map(|(x, y)| (*x.value(), *y.value())).collect();
+            bench.iter(|| {
+                for (x, y) in &i {
+                    black_box(black_box(*x) | black_box(*y));
+                }
+            });
+        });
+    }
+
+    for i in input.clone().into_iter() {
+        group.bench_with_input(format!("bitxor {:?}", &i.len()), &i, |bench, i| {
+            let i: Vec<(u32, u32)> = i.iter().map(|(x, y)| (*x.value(), *y.value())).collect();
+            bench.iter(|| {
+                for (x, y) in &i {
+                    black_box(black_box(*x) ^ black_box(*y));
+                }
+            });
+        });
+    }
+}