diff --git a/.github/workflows/rust-ci.yml b/.github/workflows/rust-ci.yml
new file mode 100644
index 00000000..63316e7b
--- /dev/null
+++ b/.github/workflows/rust-ci.yml
@@ -0,0 +1,70 @@
+name: Rust CI
+
+on: [push, pull_request]
+
+jobs:
+  build:
+    runs-on: ubuntu-latest
+
+    steps:
+      - name: Checkout code
+        uses: actions/checkout@v3
+
+      - name: Install Rust
+        uses: actions-rs/toolchain@v1
+        with:
+          toolchain: stable
+          override: true
+
+      - name: Uninstall pre-installed tools
+        run: |
+          rm -f ~/.cargo/bin/rust-analyzer
+          rm -f ~/.cargo/bin/rustfmt
+          rm -f ~/.cargo/bin/cargo-fmt
+
+      - name: Update Rust toolchain and components
+        run: |
+          rustup update
+          rustup component add rustfmt
+          rustup component add clippy
+
+      - name: Build
+        run: cargo build --verbose
+        working-directory: caledonia
+
+      - name: Verify target directory exists
+        run: |
+          echo "Checking if target directory exists and is not empty"
+          ls -la caledonia/target
+
+      - name: Run tests
+        run: cargo test --verbose
+        working-directory: caledonia
+
+#      - name: Run Clippy
+#        run: cargo clippy -- -D warnings
+#        working-directory: caledonia
+
+      - name: Check format
+        run: cargo fmt -- --check
+        working-directory: caledonia
+
+      - name: Build and test documentation
+        run: cargo doc --no-deps --verbose
+        working-directory: caledonia
+
+      - name: Cache Cargo registry
+        uses: actions/cache@v3
+        with:
+          path: ~/.cargo/registry
+          key: ${{ runner.os }}-cargo-registry-${{ hashFiles('**/Cargo.lock') }}
+          restore-keys: |
+            ${{ runner.os }}-cargo-registry-
+
+      - name: Cache Cargo build
+        uses: actions/cache@v3
+        with:
+          path: caledonia/target
+          key: ${{ runner.os }}-cargo-build-${{ hashFiles('**/Cargo.lock') }}
+          restore-keys: |
+            ${{ runner.os }}-cargo-build-
diff --git a/caledonia/.gitignore b/caledonia/.gitignore
new file mode 100644
index 00000000..9f970225
--- /dev/null
+++ b/caledonia/.gitignore
@@ -0,0 +1 @@
+target/
\ No newline at end of file
diff --git a/caledonia/Cargo.toml b/caledonia/Cargo.toml
new file mode 100644
index 00000000..5f02bb9f
--- /dev/null
+++ b/caledonia/Cargo.toml
@@ -0,0 +1,17 @@
+[package]
+name = "caledonia"
+version = "0.1.0"
+edition = "2021"
+description = "A Rust implementation of Approximate Lower Bound Arguments (ALBAs)."
+categories = ["cryptography"]
+include = ["**/*.rs", "Cargo.toml", "README.md", ".gitignore"]
+
+# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
+
+[dependencies]
+blake2 = "0.10.6"
+rand_core = "0.6.4"
+
+[dev-dependencies]
+rand = "0.8.5"
+rand_chacha = "0.3.1"
\ No newline at end of file
diff --git a/caledonia/src/bounded.rs b/caledonia/src/bounded.rs
new file mode 100644
index 00000000..9afc1293
--- /dev/null
+++ b/caledonia/src/bounded.rs
@@ -0,0 +1,566 @@
+//! Rust implementation of ALBA's bounded DFS scheme using Blake2b as hash
+//! function.
+
+extern crate core;
+use crate::utils;
+
+use std::f64::consts::E;
+
+const DATA_LENGTH: usize = 64;
+const DIGEST_SIZE: usize = 32;
+
+type Data = [u8; DATA_LENGTH];
+type Hash = [u8; DIGEST_SIZE];
+
+/// Setup input parameters
+#[derive(Debug, Clone)]
+pub struct Params {
+    /// Soundness security parameter
+    pub lambda_sec: usize,
+    /// Completeness security parameter
+    pub lambda_rel: usize,
+    /// Approximate size of set Sp to lower bound
+    pub n_p: usize,
+    /// Target lower bound
+    pub n_f: usize,
+}
+pub enum Cases {
+    /// Case where u =< λ^2
+    Small,
+    /// Case where λ^2 < u < λ^3
+    Mid,
+    /// Case where u >= λ^3
+    High,
+}
+
+impl Params {
+    /// Returns information on which case corresponds some parameter
+    pub fn which_case(&self) -> (Cases, usize) {
+        let lsec = self.lambda_sec as f64;
+        let lrel = self.lambda_rel as f64;
+        let np = self.n_p as f64;
+        let nf = self.n_f as f64;
+        let loge = E.log2();
+
+        let lognpnf = (np / nf).log2();
+        let u_f64 = (lsec + lrel.log2() + 5.0 - loge.log2()) / lognpnf;
+        let u = u_f64.ceil() as u64;
+
+        let ratio = 9.0 * np * loge / ((17 * u).pow(2) as f64);
+        let s1 = ratio - 7.0;
+        let s2 = ratio - 2.0;
+
+        if s1 < 1.0 || s2 < 1.0 {
+            return (Cases::Small, u as usize);
+        }
+
+        let lrel2 = lrel.min(s2);
+        if (u as f64) < lrel2 {
+            return (Cases::Mid, u as usize);
+        } else {
+            return (Cases::High, u as usize);
+        }
+    }
+}
+
+/// Setup output parameters
+#[derive(Debug, Clone)]
+pub struct Setup {
+    /// Approximate size of set Sp to lower bound
+    pub n_p: usize,
+    /// Proof size (in Sp elements)
+    pub u: usize,
+    /// Proof max counter
+    pub r: usize,
+    /// Proof max 2nd counter
+    pub d: usize,
+    /// Inverse of probability p_q
+    pub q: usize,
+    /// Computation bound
+    pub b: usize,
+}
+impl Setup {
+    /// Setup algorithm taking a Params as input and returning setup parameters (u,d,q)
+    pub fn new(params: &Params) -> Self {
+        let loge = E.log2();
+        fn compute_w(u: f64, l: f64) -> f64 {
+            fn factorial_check(w: f64, l: f64) -> bool {
+                let bound = 0.5f64.powf(l);
+                let factors: Vec<u64> = (1..=(w as u64 + 1)).rev().collect();
+                let mut ratio = (14.0 * w * w * (w + 2.0) * E.powf((w + 1.0) / w))
+                    / (E * (w + 2.0 - E.powf(1.0 / w)));
+
+                for f in factors {
+                    ratio /= f as f64;
+                    if ratio <= bound {
+                        return true;
+                    }
+                }
+                return false;
+            }
+            let mut w: f64 = u;
+            while !factorial_check(w, l) {
+                w += 1.0;
+            }
+            w
+        }
+
+        let n_p_f64 = params.n_p as f64;
+        let n_f_f64 = params.n_f as f64;
+        let lognpnf = (n_p_f64 / n_f_f64).log2();
+        let lambda_rel = params.lambda_rel as f64;
+        let logrel = lambda_rel.log2();
+        let lambda_sec = (params.lambda_sec as f64) + logrel;
+
+        let u_f64 = ((lambda_sec + logrel + 5.0 - loge.log2()) / lognpnf).ceil();
+        let u = u_f64 as usize;
+
+        let ratio = 9.0 * n_p_f64 * loge / ((17 * u).pow(2) as f64);
+        let s1 = ratio - 7.0;
+        let s2 = ratio - 2.0;
+
+        if s1 < 1.0 || s2 < 1.0 {
+            // Small case, ie n_p <= λ^2
+            let ln12 = (12f64).ln();
+            let d = (32.0 * ln12 * u_f64).ceil();
+            return Setup {
+                n_p: params.n_p,
+                u,
+                r: params.lambda_rel,
+                d: d as usize,
+                q: (2.0 * ln12 / d).recip().ceil() as usize,
+                b: (8.0 * (u_f64 + 1.0) * d / ln12).floor() as usize,
+            };
+        }
+        let lambda_rel2 = lambda_rel.min(s2);
+        if u_f64 < lambda_rel2 {
+            // Case 3, Theorem 14, ie  n_p >= λ^3
+            let d = (16.0 * u_f64 * (lambda_rel2 + 2.0) / loge).ceil();
+            assert!(n_p_f64 >= d * d * loge / (9.0 * (lambda_rel2 + 2.0)));
+            return Setup {
+                n_p: params.n_p,
+                u,
+                r: (lambda_rel / lambda_rel2).ceil() as usize,
+                d: d as usize,
+                q: (2.0 * (lambda_rel2 + 2.0) / (d * loge)).recip().ceil() as usize,
+                b: (((lambda_rel2 + 2.0 + u_f64.log2()) / (lambda_rel2 + 2.0))
+                    * (3.0 * u_f64 * d / 4.0)
+                    + d
+                    + u_f64)
+                    .floor() as usize,
+            };
+        } else {
+            // Case 2, Theorem 13, ie λ^3 > n_p > λ^2
+            let lambda_rel1 = lambda_rel.min(s1);
+            let lbar = (lambda_rel1 + 7.0) / loge;
+            let d = (16.0 * u_f64 * lbar).ceil();
+            assert!(n_p_f64 >= d * d / (9.0 * lbar));
+
+            let w = compute_w(u_f64, lambda_rel1);
+            return Setup {
+                n_p: params.n_p,
+                u,
+                r: (lambda_rel / lambda_rel1).ceil() as usize,
+                d: d as usize,
+                q: (2.0 * lbar / d).recip().ceil() as usize,
+                b: (((w * lbar) / d + 1.0)
+                    * E.powf(2.0 * u_f64 * w * lbar / n_p_f64 + 7.0 * u_f64 / w)
+                    * d
+                    * u_f64
+                    + d)
+                    .floor() as usize,
+            };
+        }
+    }
+}
+
+/// Round parameters
+#[derive(Debug, Clone)]
+pub struct Round {
+    /// Proof counter
+    v: usize,
+    /// Proof 2nd counter
+    t: usize,
+    // Round candidate tuple
+    s_list: Vec<Data>,
+    /// Round candidate hash
+    h: Hash,
+    /// Round candidate hash mapped to [1, n_p]
+    h_usize: usize,
+    /// Approximate size of set Sp to lower bound
+    n_p: usize,
+}
+
+impl Round {
+    /// Oracle producing a uniformly random value in [1, n_p] used for round candidates
+    /// We also return hash(data) to follow the optimization presented in Section 3.3
+    fn h1(data: Vec<Vec<u8>>, n_p: usize) -> (Hash, usize) {
+        let digest = utils::combine_hashes::<DIGEST_SIZE>(data);
+        return (digest, utils::oracle(&digest, n_p));
+    }
+
+    /// Output a round from a proof counter and n_p
+    /// Initilialises the hash with H1(t) and random value as oracle(H1(t), n_p)
+    pub fn new(v: usize, t: usize, n_p: usize) -> Round {
+        let mut data = Vec::new();
+        data.push(v.to_ne_bytes().to_vec());
+        data.push(t.to_ne_bytes().to_vec());
+        let (h, h_usize) = Round::h1(data, n_p);
+        Round {
+            v,
+            t,
+            s_list: Vec::new(),
+            h: h,
+            h_usize,
+            n_p,
+        }
+    }
+
+    /// Updates a round with an element of S_p
+    /// Replaces the hash $h$ with $h' = H1(h, s)$ and the random value as oracle(h', n_p)
+    pub fn update(r: &Round, s: Data) -> Round {
+        let mut s_list = r.s_list.clone();
+        s_list.push(s);
+        let mut data = Vec::new();
+        data.push(r.h.clone().to_vec());
+        data.push(s.to_vec());
+        let (h, h_usize) = Round::h1(data, r.n_p);
+        Round {
+            v: r.v,
+            t: r.t,
+            s_list,
+            h: h,
+            h_usize,
+            n_p: r.n_p,
+        }
+    }
+}
+
+#[derive(Debug, Clone)]
+/// Alba proof
+pub struct Proof {
+    /// Proof counter
+    r: usize,
+    /// Proof 2nd counter
+    d: usize,
+    /// Proof tuple
+    items: Vec<Data>,
+}
+
+impl Proof {
+    /// Returns a new proof
+    fn new() -> Self {
+        Proof {
+            r: 0,
+            d: 0,
+            items: Vec::new(),
+        }
+    }
+
+    /// Oracle producing a uniformly random value in [1, n_p] used for prehashing S_p
+    fn h0(setup: &Setup, v: usize, s: Data) -> usize {
+        let mut data = Vec::new();
+        data.push(v.to_ne_bytes().to_vec());
+        data.push(s.to_vec());
+        let digest = utils::combine_hashes::<DIGEST_SIZE>(data);
+        return utils::oracle(&digest, setup.n_p);
+    }
+
+    /// Oracle defined as Bernoulli(q) returning 1 with probability q and 0 otherwise
+    fn h2(setup: &Setup, r: &Round) -> bool {
+        let mut data = Vec::new();
+        data.push(r.v.to_ne_bytes().to_vec());
+        data.push(r.t.to_ne_bytes().to_vec());
+        for s in &r.s_list {
+            data.push(s.clone().to_vec());
+        }
+        let digest = utils::combine_hashes::<DIGEST_SIZE>(data);
+        return utils::oracle(&digest, setup.q) == 0;
+    }
+
+    /// Depth-first search which goes through all potential round candidates
+    /// and returns first round candidate Round{t, x_1, ..., x_u)} such that:
+    /// - for all i ∈ [0, u-1], H0(x_i+1) ∈ bins[H1(t, x_1, ..., x_i)]
+    /// - H2(t, x_0, ..., x_u) = true
+    fn dfs(
+        setup: &Setup,
+        bins: &Vec<Vec<Data>>,
+        round: &Round,
+        nb_steps: &mut usize,
+    ) -> Option<Proof> {
+        if round.s_list.len() == setup.u {
+            if Proof::h2(setup, round) {
+                let r = round.v;
+                let d = round.t;
+                let items = round.s_list.clone();
+                return Some(Proof { r, d, items });
+            } else {
+                return None;
+            }
+        }
+        let result = bins[round.h_usize].iter().find_map(|&s| {
+            if *nb_steps == setup.b {
+                return None;
+            }
+            *nb_steps += 1;
+            Self::dfs(setup, bins, &Round::update(round, s), nb_steps)
+        });
+        return result;
+    }
+
+    /// Indexed proving algorithm, returns an empty proof if no suitable
+    /// candidate is found within the setup.b steps.
+    fn prove_index(setup: &Setup, set: &Vec<Data>, v: usize) -> (usize, Option<Proof>) {
+        let mut bins: Vec<Vec<Data>> = Vec::new();
+        for _ in 1..(setup.n_p + 1) {
+            bins.push(Vec::new());
+        }
+        for &s in set.iter() {
+            bins[Proof::h0(setup, v, s)].push(s);
+        }
+        let mut nb_steps = 0;
+        for t in 1..(setup.d + 1) {
+            if nb_steps == setup.b {
+                return (0, None);
+            }
+            nb_steps += 1;
+            let round = Round::new(v, t, setup.n_p);
+            let res = Proof::dfs(setup, &bins, &round, &mut nb_steps);
+            if res.is_some() {
+                return (nb_steps, res);
+            }
+        }
+        return (nb_steps, None);
+    }
+
+    /// Alba's proving algorithm, based on a depth-first search algorithm.
+    /// Calls up to setup.r times the prove_index function and returns an empty
+    /// proof if no suitable candidate is found.
+    pub fn prove(setup: &Setup, set: &Vec<Data>) -> Self {
+        for v in 0..setup.r {
+            if let (_, Some(proof)) = Proof::prove_index(setup, set, v) {
+                return proof;
+            }
+        }
+        return Proof::new();
+    }
+
+    /// Alba's proving algorithm used for benchmarking, returning a proof as
+    /// well as the number of  steps ran to find it.
+    pub fn bench(setup: &Setup, set: &Vec<Data>) -> (usize, Self) {
+        let mut nb_steps = 0;
+        for v in 0..setup.r {
+            let (steps, opt) = Proof::prove_index(setup, set, v);
+            nb_steps += steps;
+            if let Some(proof) = opt {
+                return (nb_steps, proof);
+            }
+        }
+        return (nb_steps, Proof::new());
+    }
+
+    /// Alba's verification algorithm, follows proving algorithm by running the
+    /// same depth-first search algorithm.
+    pub fn verify(setup: &Setup, proof: Proof) -> bool {
+        if proof.d == 0 || proof.d > setup.d || proof.r > setup.r || proof.items.len() != setup.u {
+            return false;
+        }
+        let r0 = Round::new(proof.r, proof.d, setup.n_p);
+        let (b, round) = proof.items.iter().fold((true, r0), |(b, r), &s| {
+            (
+                b && r.h_usize == Proof::h0(setup, proof.r, s),
+                Round::update(&r, s),
+            )
+        });
+        return b && Proof::h2(setup, &round);
+    }
+}
+
+#[cfg(test)]
+mod tests {
+
+    use super::*;
+    use rand_chacha::ChaCha20Rng;
+    use rand_core::{RngCore, SeedableRng};
+
+    #[test]
+    fn test_params() {
+        let lambdas = [10, 80, 100, 128];
+        let pows: Vec<u32> = (2..10).collect();
+        let sps: Vec<usize> = pows.iter().map(|&i| 10_u32.pow(i) as usize).collect();
+        let ratios = [60, 66, 80, 95, 99];
+        let mut params = Vec::new();
+        for l in lambdas {
+            for &sp in &sps {
+                for r in ratios {
+                    params.push(Params {
+                        lambda_sec: l,
+                        lambda_rel: l,
+                        n_p: (sp * r) / 100,
+                        n_f: (sp * (100 - r)) / 100,
+                    })
+                }
+            }
+        }
+
+        let mut smalls = Vec::new();
+        let mut mids = Vec::new();
+        let mut highs = Vec::new();
+        for p in params {
+            match Params::which_case(&p) {
+                (Cases::Small, u) => smalls.push((p.clone(), u)),
+                (Cases::Mid, u) => mids.push((p.clone(), u)),
+                (Cases::High, u) => highs.push((p.clone(), u)),
+            }
+        }
+
+        println!("------------ Small cases");
+        for s in smalls {
+            println!("{:?}", s);
+        }
+        println!("\n------------ Mid cases");
+        for s in mids {
+            println!("{:?}", s);
+        }
+        println!("\n------------ High cases");
+        for s in highs {
+            println!("{:?}", s);
+        }
+    }
+
+    #[test]
+    fn test_verify() {
+        let mut rng = ChaCha20Rng::from_seed(Default::default());
+        let nb_tests = 1_000;
+        let set_size = 1_000;
+        for _t in 0..nb_tests {
+            let seed = rng.next_u32().to_ne_bytes().to_vec();
+            let s_p = utils::gen_items::<DATA_LENGTH>(seed, set_size);
+            let params = Params {
+                lambda_sec: 10,
+                lambda_rel: 10,
+                n_p: 80,
+                n_f: 20,
+            };
+            let setup = Setup::new(&params);
+            let proof = Proof::prove(&setup, &s_p);
+            assert!(Proof::verify(&setup, proof.clone()));
+            let proof_0 = Proof {
+                r: proof.r,
+                d: 0,
+                items: proof.items.clone(),
+            };
+            assert!(!Proof::verify(&setup, proof_0));
+            let proof_d = Proof {
+                r: proof.r,
+                d: proof.d.wrapping_add(1),
+                items: proof.items.clone(),
+            };
+            assert!(!Proof::verify(&setup, proof_d));
+            let proof_r = Proof {
+                r: proof.r.wrapping_add(1),
+                d: proof.d,
+                items: proof.items.clone(),
+            };
+            assert!(!Proof::verify(&setup, proof_r));
+            let proof_item = Proof {
+                r: proof.r,
+                d: proof.d,
+                items: Vec::new(),
+            };
+            assert!(!Proof::verify(&setup, proof_item));
+            let mut wrong_items = proof.items.clone();
+            let last_item = wrong_items.pop().unwrap();
+            let mut penultimate_item = wrong_items.pop().unwrap();
+            let proof_itembis = Proof {
+                r: proof.r,
+                d: proof.d,
+                items: wrong_items.clone(),
+            };
+            assert!(!Proof::verify(&setup, proof_itembis));
+            // Modifying the penultimate item to check correctness of H1 check and not H2
+            penultimate_item[0] = penultimate_item[0].wrapping_add(42u8);
+            wrong_items.push(penultimate_item);
+            wrong_items.push(last_item);
+            let proof_itembis = Proof {
+                r: proof.r,
+                d: proof.d,
+                items: wrong_items.clone(),
+            };
+            assert!(!Proof::verify(&setup, proof_itembis));
+        }
+    }
+
+    #[test]
+    fn test_prove() {
+        use std::time::Instant;
+        let npnf = [
+            // (99_000, 1_000),  // high
+            // (95_000, 5_000),  // medium
+            // (80_000, 20_000), // low
+            // (66_000, 34_000),
+            (60_000, 40_000),
+        ];
+        let lambdas = [80];
+        let nb_tests = 100;
+        let mut rng = ChaCha20Rng::from_seed([0u8; 32]);
+        for (n_p, n_f) in npnf {
+            for lambda in lambdas {
+                let mut u = 0;
+                let mut time_setup = 0;
+                let mut time_prove = 0;
+                let mut time_verify = 0;
+                let mut max_bench: usize = 0;
+                let mut mean_bench: usize = 0;
+                let mut max_retrial: usize = 0;
+                let mut mean_retrial: usize = 0;
+                for _t in 0..nb_tests {
+                    let seed_u32 = rng.next_u32();
+                    let seed = seed_u32.to_ne_bytes().to_vec();
+                    let s_p: Vec<Data> = utils::gen_items::<DATA_LENGTH>(seed, n_p);
+                    let params = Params {
+                        lambda_sec: lambda,
+                        lambda_rel: lambda,
+                        n_p,
+                        n_f,
+                    };
+                    // Setup
+                    let start_setup = Instant::now();
+                    let setup = Setup::new(&params);
+                    let end_setup = start_setup.elapsed();
+                    time_setup += end_setup.as_nanos();
+                    u = setup.u;
+                    // Prove
+                    let start_prove = Instant::now();
+                    let (steps, proof) = Proof::bench(&setup, &s_p);
+                    let end_prove = start_prove.elapsed();
+                    time_prove += end_prove.as_nanos();
+                    max_bench = std::cmp::max(max_bench, steps);
+                    mean_bench += steps;
+                    max_retrial = std::cmp::max(max_retrial, proof.r);
+                    mean_retrial += proof.r;
+                    // Verify
+                    let start_verify = Instant::now();
+                    let b = Proof::verify(&setup, proof.clone());
+                    let end_verify = start_verify.elapsed();
+                    time_verify += end_verify.as_nanos();
+                    assert!(b);
+                }
+                println!(
+                    "(n_p={}, n_f={}, λ={}): \t u={}, \t setup:{}, \t prove:{}, \t verify:{}, \t max steps:{}, \t mean steps:{}, \t max retrial:{}, \t mean retrial:{}",
+                    utils::format_nb(n_p),
+                    utils::format_nb(n_f),
+                    utils::format_nb(lambda),
+                    utils::format_nb(u),
+                    utils::format_time(time_setup / nb_tests),
+                    utils::format_time(time_prove / nb_tests),
+                    utils::format_time(time_verify / nb_tests),
+                    utils::format_nb(max_bench),
+                    utils::format_nb((mean_bench as u128 / nb_tests) as usize),
+                    max_retrial,
+                    mean_retrial as u128 / nb_tests,
+                );
+            }
+        }
+    }
+}
diff --git a/caledonia/src/lib.rs b/caledonia/src/lib.rs
new file mode 100644
index 00000000..597b3ca7
--- /dev/null
+++ b/caledonia/src/lib.rs
@@ -0,0 +1,17 @@
+//! Approximate Lower Bound Arguments (ALBA, <https://eprint.iacr.org/2023/1655.pdf>)
+//!
+//! Alba is a generic protocol to prove succinctly a lower bound of the size of
+//! a, potentially weighted, set. Say we have a set Sp of size |Sp| >= $n_p$,
+//! and a lower bound $n_f$ < $n_p$ of it we want to prove. Alba gives us a
+//! method to generate a proof of knowledge of this bound by finding the
+//! smallest subset of Sp of size $u$ to convince a verifier.
+//! The paper presents several schemes and optimizations. The basic scheme is
+//! enhanced in the "prehashed" version thanks to sorting Sp with a balls and
+//! bins sorting algorithm reducing the number of hashes done per round. A
+//! lottery scheme is also introduced to support Alba in a decentralised
+//! settings as well as a modification to use PRF in the CRS settings instead of
+//! using the ROM.
+pub mod utils;
+
+pub mod bounded;
+pub mod prehashed;
diff --git a/caledonia/src/prehashed.rs b/caledonia/src/prehashed.rs
new file mode 100644
index 00000000..9898bc26
--- /dev/null
+++ b/caledonia/src/prehashed.rs
@@ -0,0 +1,331 @@
+//! Rust implementation of ALBA's prehashed scheme using Blake2b as hash
+//! function, working for big sets.
+
+extern crate core;
+use crate::utils;
+
+use std::f64::consts::E;
+
+const DATA_LENGTH: usize = 32;
+const DIGEST_SIZE: usize = 32;
+
+type Data = [u8; DATA_LENGTH];
+type Hash = [u8; DIGEST_SIZE];
+
+/// Setup input parameters
+#[derive(Debug, Clone)]
+pub struct Params {
+    /// Soundness security parameter
+    pub lambda_sec: usize,
+    /// Completeness security parameter
+    pub lambda_rel: usize,
+    /// Approximate size of set Sp to lower bound
+    pub n_p: usize,
+    /// Target lower bound
+    pub n_f: usize,
+}
+/// Setup output parameters
+#[derive(Debug, Clone)]
+pub struct Setup {
+    /// Approximate size of set Sp to lower bound
+    pub n_p: usize,
+    /// Proof size (in Sp elements)
+    pub u: usize,
+    /// Proof max counter
+    pub d: usize,
+    /// Inverse of probability p_q
+    pub q: usize,
+}
+impl Setup {
+    /// Setup algorithm taking a Params as input and returning setup parameters (u,d,q)
+    pub fn new(params: &Params) -> Self {
+        let e = E;
+        let log_2 = |x: f64| x.log2();
+        let loge = log_2(e);
+        let logloge = log_2(loge);
+        let log3 = log_2(3.0);
+
+        let n_p_f64 = params.n_p as f64;
+        let n_f_f64 = params.n_f as f64;
+        let lognpnf = log_2(n_p_f64 / n_f_f64);
+        let lambda_sec = params.lambda_sec as f64;
+        let lambda_rel = params.lambda_rel as f64;
+
+        // We define the parameters according to Section 3.2, Corrollary 2
+        let u_f64 = (lambda_sec + log_2(lambda_rel + log3) + 1.0 - logloge) / lognpnf;
+        let u = u_f64.ceil() as usize;
+        let d_f64 = 16.0 * u_f64 * (lambda_rel + log3) / loge;
+        let d = d_f64.ceil() as usize;
+        let q = (2.0 * (lambda_rel + log3) / (d_f64 * loge)).recip().ceil() as usize;
+
+        let check = ((d_f64 * d_f64 * loge) / (9.0 * (lambda_rel + log3))).ceil() as usize;
+        assert!(params.n_p >= check);
+
+        Setup {
+            n_p: params.n_p,
+            u,
+            d,
+            q,
+        }
+    }
+}
+
+/// Round parameters
+#[derive(Debug, Clone)]
+pub struct Round {
+    /// Proof counter
+    t: usize,
+    // Round candidate tuple
+    s_list: Vec<Data>,
+    /// Round candidate hash
+    h: Vec<u8>,
+    /// Round candidate hash mapped to [1, n_p]
+    h_usize: usize,
+    /// Approximate size of set Sp to lower bound
+    n_p: usize,
+}
+
+impl Round {
+    /// Oracle producing a uniformly random value in [1, n_p] used for round candidates
+    /// We also return hash(data) to follow the optimization presented in Section 3.3
+    fn h1(data: Vec<Vec<u8>>, n_p: usize) -> (Hash, usize) {
+        let digest = utils::combine_hashes(data);
+        return (digest, utils::oracle(&digest, n_p));
+    }
+
+    /// Output a round from a proof counter and n_p
+    /// Initilialises the hash with H1(t) and random value as oracle(H1(t), n_p)
+    pub fn new(t: usize, n_p: usize) -> Round {
+        let data = [t.to_ne_bytes().to_vec()].to_vec();
+        let (h, h_usize) = Round::h1(data, n_p);
+        Round {
+            t,
+            s_list: Vec::new(),
+            h: h.to_vec(),
+            h_usize,
+            n_p,
+        }
+    }
+
+    /// Updates a round with an element of S_p
+    /// Replaces the hash $h$ with $h' = H1(h, s)$ and the random value as oracle(h', n_p)
+    pub fn update(r: &Round, s: Data) -> Round {
+        let mut s_list = r.s_list.clone();
+        s_list.push(s);
+        let mut data = Vec::new();
+        data.push(r.h.clone());
+        data.push(s.to_vec());
+        let (h, h_usize) = Round::h1(data, r.n_p);
+        Round {
+            t: r.t,
+            s_list,
+            h: h.to_vec(),
+            h_usize,
+            n_p: r.n_p,
+        }
+    }
+}
+
+#[derive(Debug, Clone)]
+/// Alba proof
+pub struct Proof {
+    /// Proof counter
+    d: usize,
+    /// Proof tuple
+    items: Vec<Data>,
+}
+
+impl Proof {
+    /// Oracle producing a uniformly random value in [1, n_p] used for prehashing S_p
+    // TODO: We also return hash(data) to follow the optimization presented in Section 3.3
+    fn h0(setup: &Setup, s: Data) -> usize {
+        let mut data = Vec::new();
+        data.push(s.to_vec());
+        let digest = utils::combine_hashes::<DIGEST_SIZE>(data);
+        return utils::oracle(&digest, setup.n_p);
+    }
+
+    /// Oracle defined as Bernoulli(q) returning 1 with probability q and 0 otherwise
+    fn h2(setup: &Setup, r: &Round) -> bool {
+        let mut data = Vec::new();
+        data.push(r.t.to_ne_bytes().to_vec());
+        for s in &r.s_list {
+            data.push(s.clone().to_vec());
+        }
+        let digest = utils::combine_hashes::<DIGEST_SIZE>(data);
+        return utils::oracle(&digest, setup.q) == 0;
+    }
+
+    /// Depth-first search which goes through all potential round candidates
+    /// and returns first round candidate Round{t, x_1, ..., x_u)} such that:
+    /// - for all i ∈ [0, u-1], H0(x_i+1) ∈ bins[H1(t, x_1, ..., x_i)]
+    /// - H2(t, x_0, ..., x_u) = true
+    fn dfs(setup: &Setup, bins: &Vec<Vec<Data>>, round: &Round) -> Option<Proof> {
+        if round.s_list.len() == setup.u {
+            if Proof::h2(setup, round) {
+                let d = round.t;
+                let items = round.s_list.clone();
+                return Some(Proof { d, items });
+            } else {
+                return None;
+            }
+        }
+        let result = bins[round.h_usize]
+            .iter()
+            .find_map(|&s| Self::dfs(setup, bins, &Round::update(round, s)));
+        return result;
+    }
+
+    /// Alba's proving algorithm, based on a depth-first search algorithm.
+    /// Returns an empty proof if no suitable candidate is found.
+    pub fn prove(setup: &Setup, set: &Vec<Data>) -> Self {
+        let mut bins: Vec<Vec<Data>> = Vec::new();
+        for _ in 1..(setup.n_p + 1) {
+            bins.push(Vec::new());
+        }
+        for &s in set.iter() {
+            bins[Proof::h0(setup, s)].push(s);
+        }
+
+        for t in 1..(setup.d + 1) {
+            let round = Round::new(t, setup.n_p);
+            if let Some(proof) = Proof::dfs(setup, &bins, &round) {
+                return proof;
+            };
+        }
+
+        return Proof {
+            d: 0,
+            items: Vec::new(),
+        };
+    }
+
+    /// Alba's verification algorithm, follows proving algorithm by running the
+    /// same depth-first search algorithm.
+    pub fn verify(setup: &Setup, proof: Proof) -> bool {
+        if proof.d == 0 || proof.d > setup.d || proof.items.len() != setup.u {
+            return false;
+        }
+        let r0 = Round::new(proof.d, setup.n_p);
+        let (b, round) = proof.items.iter().fold((true, r0), |(b, r), &s| {
+            (b && r.h_usize == Proof::h0(setup, s), Round::update(&r, s))
+        });
+        return b && Proof::h2(setup, &round);
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use rand_chacha::ChaCha20Rng;
+    use rand_core::{RngCore, SeedableRng};
+
+    #[test]
+    fn test_verify() {
+        let mut rng = ChaCha20Rng::from_seed(Default::default());
+        let nb_tests = 1_000;
+        let set_size = 1_000;
+        for _t in 0..nb_tests {
+            let seed = rng.next_u32().to_ne_bytes().to_vec();
+            let s_p = utils::gen_items::<DATA_LENGTH>(seed, set_size);
+            let params = Params {
+                lambda_sec: 10,
+                lambda_rel: 10,
+                n_p: 800,
+                n_f: 2,
+            };
+            let setup = Setup::new(&params);
+            let proof = Proof::prove(&setup, &s_p);
+            assert!(Proof::verify(&setup, proof.clone()));
+            let proof_0 = Proof {
+                d: 0,
+                items: proof.items.clone(),
+            };
+            assert!(!Proof::verify(&setup, proof_0));
+            let proof_d = Proof {
+                d: proof.d.wrapping_add(1),
+                items: proof.items.clone(),
+            };
+            assert!(!Proof::verify(&setup, proof_d));
+            let proof_item = Proof {
+                d: proof.d,
+                items: Vec::new(),
+            };
+            assert!(!Proof::verify(&setup, proof_item));
+            let mut wrong_items = proof.items.clone();
+            let last_item = wrong_items.pop().unwrap();
+            let mut penultimate_item = wrong_items.pop().unwrap();
+            let proof_itembis = Proof {
+                d: proof.d,
+                items: wrong_items.clone(),
+            };
+            assert!(!Proof::verify(&setup, proof_itembis));
+            // Modifying the penultimate item to check correctness of H1 check and not H2
+            penultimate_item[0] = penultimate_item[0].wrapping_add(1);
+            wrong_items.push(penultimate_item);
+            wrong_items.push(last_item);
+            let proof_itembis = Proof {
+                d: proof.d,
+                items: wrong_items,
+            };
+            assert!(!Proof::verify(&setup, proof_itembis));
+        }
+    }
+
+    #[test]
+    fn test_prove() {
+        use std::time::Instant;
+        let npnf = [(1_000, 8)];
+        // Other working tests: (10_000, 2_000), (100_000, 60_000)
+        let lambdas = [5, 10];
+        let nb_tests = 100;
+        let mut rng = ChaCha20Rng::from_seed([0u8; 32]);
+        for (n_p, n_f) in npnf {
+            let mut u = 0;
+            let mut time_setup = 0;
+            let mut time_prove = 0;
+            let mut time_verify = 0;
+            for lambda in lambdas {
+                for _t in 0..nb_tests {
+                    let seed_u32 = rng.next_u32();
+                    let seed = seed_u32.to_ne_bytes().to_vec();
+                    let s_p: Vec<Data> = utils::gen_items::<DATA_LENGTH>(seed, n_p);
+                    let params = Params {
+                        lambda_sec: lambda,
+                        lambda_rel: lambda,
+                        n_p,
+                        n_f,
+                    };
+                    // Setup
+                    let start_setup = Instant::now();
+                    let setup = Setup::new(&params);
+                    let end_setup = start_setup.elapsed();
+                    time_setup += end_setup.as_nanos();
+                    u = setup.u;
+                    // Prove
+                    let start_prove = Instant::now();
+                    let proof = Proof::prove(&setup, &s_p);
+                    let end_prove = start_prove.elapsed();
+                    time_prove += end_prove.as_nanos();
+                    // Verify
+                    let start_verify = Instant::now();
+                    let b = Proof::verify(&setup, proof.clone());
+                    let end_verify = start_verify.elapsed();
+                    time_verify += end_verify.as_nanos();
+                    assert!(b);
+                }
+
+                println!(
+                    "(n_p={}, n_f={}, λ={}): \t u={}, \t setup:{}, \t prove:{}, \t verify:{}",
+                    utils::format_nb(n_p),
+                    utils::format_nb(n_f),
+                    lambda,
+                    utils::format_nb(u),
+                    utils::format_time(time_setup / nb_tests),
+                    utils::format_time(time_prove / nb_tests),
+                    utils::format_time(time_verify / nb_tests)
+                );
+            }
+        }
+    }
+}
diff --git a/caledonia/src/utils.rs b/caledonia/src/utils.rs
new file mode 100644
index 00000000..13cc4799
--- /dev/null
+++ b/caledonia/src/utils.rs
@@ -0,0 +1,123 @@
+use blake2::digest::{Update, VariableOutput};
+use blake2::Blake2bVar;
+use std::cmp::min;
+
+// Helper functions
+fn mod_non_power_of_2(hash: &[u8], n: usize) -> usize {
+    let epsilon_fail: usize = 1 << 40; // roughly 1 in 10 billion
+    let k: usize = log_base2(n * epsilon_fail);
+    let k_prime: usize = 1 << k;
+    let d: usize = k_prime.div_ceil(n);
+
+    let i = mod_power_of_2(hash, k_prime);
+
+    if i >= d * n {
+        panic!("failed: i = {}, d = {}, n = {}, k = {}", i, d, n, k);
+    } else {
+        i % n
+    }
+}
+
+fn mod_power_of_2(hash: &[u8], n: usize) -> usize {
+    let r = from_bytes_le(hash);
+    (n - 1) & r
+}
+
+fn log_base2(x: usize) -> usize {
+    usize::BITS as usize - x.leading_zeros() as usize - 1
+}
+
+fn from_bytes_le(bytes: &[u8]) -> usize {
+    let mut array = [0u8; 8];
+    let bytes = &bytes[..min(8, bytes.len())];
+    array[..bytes.len()].copy_from_slice(bytes);
+    usize::from_le_bytes(array)
+}
+
+/// Return a 32-byte hash of the given data
+pub fn hash_bytes<const N: usize>(data: &[u8]) -> [u8; N] {
+    let mut hasher = Blake2bVar::new(N).expect("Failed to construct hasher!");
+    hasher.update(data);
+    let mut buf = [0u8; N];
+    hasher
+        .finalize_variable(&mut buf)
+        .expect("Failed to finalize hashing");
+    buf
+}
+
+/// Return 32-byte hash of the given list of data
+pub fn combine_hashes<const N: usize>(hash_list: Vec<Vec<u8>>) -> [u8; N] {
+    let mut hasher = Blake2bVar::new(N).expect("Failed to construct hasher!");
+    for data in hash_list.iter() {
+        hasher.update(data);
+    }
+    let mut buf = [0u8; N];
+    hasher
+        .finalize_variable(&mut buf)
+        .expect("Failed to finalize hashing");
+    buf
+}
+
+pub fn oracle(hash: &[u8], n: usize) -> usize {
+    if n.is_power_of_two() {
+        mod_power_of_2(hash, n)
+    } else {
+        mod_non_power_of_2(hash, n)
+    }
+}
+
+/// Generate a set of items for given set size
+/// Items are generated by hashing the current index
+pub fn gen_items<const N: usize>(seed: Vec<u8>, set_size: usize) -> Vec<[u8; N]> {
+    let mut s_p = Vec::with_capacity(set_size);
+    for b in 0..set_size {
+        let mut data = Vec::new();
+        data.push(seed.clone());
+        data.push(b.to_ne_bytes().to_vec());
+        let item = combine_hashes::<N>(data);
+        s_p.push(item);
+    }
+    s_p
+}
+
+pub fn format_time(nanos: u128) -> String {
+    let mut time = nanos;
+    let bounds = [1000, 1000, 1000, 60, 60, 60];
+    let units = ["ns", "μs", "ms", "s", "min", "h"];
+    for (&bound, &unit) in bounds.iter().zip(units.iter()) {
+        if time < bound {
+            return time.to_string() + unit;
+        }
+        time = time / bound;
+    }
+    (time * 60).to_string() + "h"
+}
+
+pub fn format_nb(x: usize) -> String {
+    let mut y = x;
+    let mut s = String::new();
+    let mut b = true;
+    while y / 1000 != 0 {
+        let to_add = (y % 1000).to_string();
+        let preppend = "0".repeat(3 - to_add.len()) + &to_add;
+        let append = if b { "" } else { &("_".to_string() + &s) };
+        s = preppend + append;
+        b = false;
+        y = y / 1000;
+    }
+    if b {
+        y.to_string()
+    } else {
+        (y % 1000).to_string() + "_" + &s
+    }
+}
+
+#[cfg(test)]
+mod tests {
+
+    #[test]
+    fn test_oracle() {
+        // Test distribution of oracle
+        assert!(true);
+    }
+}