diff --git a/examples/bls-multisig/python/.gitignore b/examples/bls-multisig/python/.gitignore
new file mode 100644
index 0000000..21d0b89
--- /dev/null
+++ b/examples/bls-multisig/python/.gitignore
@@ -0,0 +1 @@
+.venv/
diff --git a/examples/bls-multisig/python/README.md b/examples/bls-multisig/python/README.md
new file mode 100644
index 0000000..bd10135
--- /dev/null
+++ b/examples/bls-multisig/python/README.md
@@ -0,0 +1,19 @@
+# Python Alphanet Multisig
+
+This example demonstrates an integration of [BlsMultisig](../../../src/BLSMultisig.sol) with Python.
+
+## Running the example
+
+To run the example, you will need to install the required dependencies:
+
+```shell
+pip install web3 py_ecc
+```
+
+Then, you can run the example by executing the following command:
+
+```shell
+python multisig.py
+```
+
+This will spin up an Anvil instance in Alphanet mode, deploy the multisig contract and execute a simple operation signed by random BLS keys.
\ No newline at end of file
diff --git a/examples/bls-multisig/python/multisig.py b/examples/bls-multisig/python/multisig.py
new file mode 100644
index 0000000..6399437
--- /dev/null
+++ b/examples/bls-multisig/python/multisig.py
@@ -0,0 +1,105 @@
+from web3 import AsyncWeb3
+import pathlib
+import asyncio
+import json
+import subprocess
+import random
+from py_ecc.bls import G2Basic
+from py_ecc.bls import g2_primitives
+import eth_abi
+
+Fp = tuple[int, int]
+Fp2 = tuple[Fp, Fp]
+G1Point = tuple[Fp, Fp]
+G2Point = tuple[Fp2, Fp2]
+Operation = tuple[str, str, int, int]
+
+
+def fp_from_int(x: int) -> Fp:
+    b = x.to_bytes(64, "big")
+    return (int.from_bytes(b[:32], "big"), int.from_bytes(b[32:], "big"))
+
+
+def generate_keys(num: int) -> list[tuple[G1Point, int]]:
+    keypairs = []
+    for _ in range(num):
+        sk = random.randint(0, 10**30)
+        pk_point = g2_primitives.pubkey_to_G1(G2Basic.SkToPk(sk))
+
+        pk = (fp_from_int(int(pk_point[0])), fp_from_int(int(pk_point[1])))
+
+        keypairs.append((pk, sk))
+
+    keypairs.sort()
+
+    return keypairs
+
+
+def sign_operation(sks: list[int], operation: Operation) -> G2Point:
+    encoded = eth_abi.encode(["(address,bytes,uint256,uint256)"], [operation])
+
+    signatures = []
+    for sk in sks:
+        signatures.append(G2Basic.Sign(sk, encoded))
+
+    aggregated = g2_primitives.signature_to_G2(G2Basic.Aggregate(signatures))
+
+    signature = (
+        (fp_from_int(aggregated[0].coeffs[0]), fp_from_int(aggregated[0].coeffs[1])),
+        (fp_from_int(aggregated[1].coeffs[0]), fp_from_int(aggregated[1].coeffs[1])),
+    )
+
+    return signature
+
+
+async def main():
+    bls_multisig_artifact = json.load(
+        open(pathlib.Path(__file__).parent.parent.parent.parent / "out/BLSMultisig.sol/BLSMultisig.json")
+    )
+
+    web3 = AsyncWeb3(AsyncWeb3.AsyncHTTPProvider("http://localhost:8545"))
+
+    bytecode = bls_multisig_artifact["bytecode"]["object"]
+    abi = bls_multisig_artifact["abi"]
+    BlsMultisig = web3.eth.contract(abi=abi, bytecode=bytecode)
+
+    signer = (await web3.eth.accounts)[0]
+
+    # generate 100 BLS keys
+    keypairs = generate_keys(100)
+    pks = list(map(lambda x: x[0], keypairs))
+
+    # deploy the multisig contract with generated signers and threshold of 50
+    tx = await BlsMultisig.constructor(pks, 50).transact({"from": signer})
+    receipt = await web3.eth.wait_for_transaction_receipt(tx)
+    multisig = BlsMultisig(receipt.contractAddress)
+
+    # fund the multisig
+    hash = await web3.eth.send_transaction({"from": signer, "to": multisig.address, "value": 10**18})
+    await web3.eth.wait_for_transaction_receipt(hash)
+
+    # create an operation transferring 1 eth to zero address
+    operation: Operation = ("0x0000000000000000000000000000000000000000", bytes(), 10**18, 0)
+
+    # choose 50 random signers that will sign the operation
+    signers_subset = sorted(random.sample(keypairs, 50))
+
+    pks = list(map(lambda x: x[0], signers_subset))
+    sks = list(map(lambda x: x[1], signers_subset))
+
+    # create aggregated signature for operation
+    signature = sign_operation(sks, operation)
+
+    # execute the operation
+    tx = await multisig.functions.verifyAndExecute((operation, pks, signature)).transact({"from": signer})
+    receipt = await web3.eth.wait_for_transaction_receipt(tx)
+
+    assert receipt.status == 1
+
+
+if __name__ == "__main__":
+    try:
+        anvil = subprocess.Popen(["anvil", "--alphanet"], stdout=subprocess.PIPE)
+        asyncio.run(main())
+    finally:
+        anvil.terminate()
diff --git a/examples/bls-multisig/rust/.gitignore b/examples/bls-multisig/rust/.gitignore
new file mode 100644
index 0000000..2c96eb1
--- /dev/null
+++ b/examples/bls-multisig/rust/.gitignore
@@ -0,0 +1,2 @@
+target/
+Cargo.lock
diff --git a/examples/bls-multisig/rust/Cargo.toml b/examples/bls-multisig/rust/Cargo.toml
new file mode 100644
index 0000000..17f8510
--- /dev/null
+++ b/examples/bls-multisig/rust/Cargo.toml
@@ -0,0 +1,15 @@
+[package]
+name = "alphanet-bls-multisig"
+version = "0.1.0"
+edition = "2021"
+
+[dependencies]
+alloy = { version = "0.4", features = [
+    "providers",
+    "contract",
+    "sol-types",
+    "node-bindings",
+] }
+tokio = { version = "1", features = ["full"] }
+blst = "0.3"
+rand = "0.8"
diff --git a/examples/bls-multisig/rust/src/main.rs b/examples/bls-multisig/rust/src/main.rs
new file mode 100644
index 0000000..17787a1
--- /dev/null
+++ b/examples/bls-multisig/rust/src/main.rs
@@ -0,0 +1,100 @@
+use alloy::{primitives::U256, providers::ProviderBuilder, sol, sol_types::SolValue};
+use blst::min_pk::{AggregateSignature, SecretKey, Signature};
+use rand::RngCore;
+use BLS::G2Point;
+
+sol! {
+    #[derive(Debug, Default, PartialEq, Eq, PartialOrd, Ord)]
+    #[sol(rpc)]
+    BLSMultisig,
+    "../../../out/BLSMultisig.sol/BLSMultisig.json"
+}
+
+impl From<[u8; 96]> for BLS::G1Point {
+    fn from(value: [u8; 96]) -> Self {
+        let mut data = [0u8; 128];
+        data[16..64].copy_from_slice(&value[0..48]);
+        data[80..128].copy_from_slice(&value[48..96]);
+
+        BLS::G1Point::abi_decode(&data, false).unwrap()
+    }
+}
+
+impl From<[u8; 192]> for BLS::G2Point {
+    fn from(value: [u8; 192]) -> Self {
+        let mut data = [0u8; 256];
+        data[16..64].copy_from_slice(&value[48..96]);
+        data[80..128].copy_from_slice(&value[0..48]);
+        data[144..192].copy_from_slice(&value[144..192]);
+        data[208..256].copy_from_slice(&value[96..144]);
+
+        BLS::G2Point::abi_decode(&data, false).unwrap()
+    }
+}
+
+/// Generates `num` BLS keys and returns them as a tuple of secret keys and public keys, sorted by public key.
+fn generate_keys(num: usize) -> (Vec<SecretKey>, Vec<BLS::G1Point>) {
+    let mut rng = rand::thread_rng();
+    let mut keys = Vec::with_capacity(num);
+
+    for _ in 0..num {
+        let mut ikm = [0u8; 32];
+        rng.fill_bytes(&mut ikm);
+
+        let sk = SecretKey::key_gen(&ikm, &[]).unwrap();
+        let pk: BLS::G1Point = sk.sk_to_pk().serialize().into();
+
+        keys.push((sk, pk));
+    }
+
+    keys.sort_by(|(_, pk1), (_, pk2)| pk1.cmp(pk2));
+
+    keys.into_iter().unzip()
+}
+
+/// Signs a message with the provided keys and returns the aggregated signature.
+fn sign_message(keys: &[SecretKey], msg: &[u8]) -> G2Point {
+    let mut sigs = Vec::new();
+
+    // create individual signatures
+    for key in keys {
+        let sig = key.sign(msg, b"BLS_SIG_BLS12381G2_XMD:SHA-256_SSWU_RO_NUL_", &[]);
+        sigs.push(sig);
+    }
+
+    let agg_sig = Signature::from_aggregate(
+        &AggregateSignature::aggregate(sigs.iter().collect::<Vec<_>>().as_slice(), false).unwrap(),
+    );
+
+    agg_sig.serialize().into()
+}
+
+#[tokio::main]
+pub async fn main() {
+    let provider = ProviderBuilder::new().on_anvil_with_config(|config| config.arg("--alphanet"));
+
+    let (keys, signers) = generate_keys(100);
+
+    let multisig = BLSMultisig::deploy(provider, signers.clone(), U256::from(1))
+        .await
+        .unwrap();
+
+    let operation = BLSMultisig::Operation::default();
+
+    let signature = sign_message(&keys, &operation.abi_encode());
+
+    let receipt = multisig
+        .verifyAndExecute(BLSMultisig::SignedOperation {
+            operation: operation.clone(),
+            signers,
+            signature,
+        })
+        .send()
+        .await
+        .unwrap()
+        .get_receipt()
+        .await
+        .unwrap();
+
+    assert!(receipt.status());
+}
diff --git a/src/BLSMultisig.sol b/src/BLSMultisig.sol
new file mode 100644
index 0000000..2b37284
--- /dev/null
+++ b/src/BLSMultisig.sol
@@ -0,0 +1,112 @@
+// SPDX-License-Identifier: MIT
+pragma solidity ^0.8.23;
+
+import {BLS} from "./sign/BLS.sol";
+
+/// @notice BLS-powered multisignature wallet, demonstrating the use of
+/// aggregated BLS signatures for verification
+/// @dev This is for demonstration purposes only, do not use in production. This contract does
+/// not include protection from rogue public-key attacks.
+contract BLSMultisig {
+    /// @notice Public keys of signers. This may contain a pre-aggregated
+    /// public keys for common sets of signers as well.
+    mapping(bytes32 => bool) public signers;
+
+    struct Operation {
+        address to;
+        bytes data;
+        uint256 value;
+        uint256 nonce;
+    }
+
+    struct SignedOperation {
+        Operation operation;
+        BLS.G1Point[] signers;
+        BLS.G2Point signature;
+    }
+
+    /// @notice The negated generator point in G1 (-P1). Used during pairing as a first G1 point.
+    BLS.G1Point NEGATED_G1_GENERATOR = BLS.G1Point(
+        BLS.Fp(
+            31827880280837800241567138048534752271,
+            88385725958748408079899006800036250932223001591707578097800747617502997169851
+        ),
+        BLS.Fp(
+            22997279242622214937712647648895181298,
+            46816884707101390882112958134453447585552332943769894357249934112654335001290
+        )
+    );
+
+    /// @notice The number of signatures required to execute an operation.
+    uint256 public threshold;
+
+    /// @notice Nonce used for replay protection.
+    uint256 public nonce;
+
+    constructor(BLS.G1Point[] memory _signers, uint256 _threshold) {
+        for (uint256 i = 0; i < _signers.length; i++) {
+            signers[keccak256(abi.encode(_signers[i]))] = true;
+        }
+        threshold = _threshold;
+    }
+
+    /// @notice Maps an operation to a point on G2 which needs to be signed.
+    function getOperationPoint(Operation memory op) public view returns (BLS.G2Point memory) {
+        return BLS.hashToCurveG2(abi.encode(op));
+    }
+
+    /// @notice Accepts an operation signed by a subset of the signers and executes it
+    function verifyAndExecute(SignedOperation memory operation) public {
+        require(operation.operation.nonce == nonce++, "invalid nonce");
+        require(operation.signers.length >= threshold, "not enough signers");
+
+        BLS.G1Point memory aggregatedSigner;
+
+        for (uint256 i = 0; i < operation.signers.length; i++) {
+            BLS.G1Point memory signer = operation.signers[i];
+            require(signers[keccak256(abi.encode(signer))], "invalid signer");
+
+            if (i == 0) {
+                aggregatedSigner = signer;
+            } else {
+                aggregatedSigner = BLS.G1Add(aggregatedSigner, signer);
+                require(_comparePoints(operation.signers[i - 1], signer), "signers not sorted");
+            }
+        }
+
+        BLS.G1Point[] memory g1Points = new BLS.G1Point[](2);
+        BLS.G2Point[] memory g2Points = new BLS.G2Point[](2);
+
+        g1Points[0] = NEGATED_G1_GENERATOR;
+        g1Points[1] = aggregatedSigner;
+
+        g2Points[0] = operation.signature;
+        g2Points[1] = getOperationPoint(operation.operation);
+
+        // verify signature
+        require(BLS.Pairing(g1Points, g2Points), "invalid signature");
+
+        // execute operation
+        Operation memory op = operation.operation;
+        (bool success,) = op.to.call{value: op.value}(op.data);
+        require(success, "execution failed");
+    }
+
+    /// @notice Returns true if X coordinate of the first point is lower than the X coordinate of the second point.
+    function _comparePoints(BLS.G1Point memory a, BLS.G1Point memory b) internal pure returns (bool) {
+        BLS.Fp memory aX = a.x;
+        BLS.Fp memory bX = b.x;
+
+        if (aX.a < bX.a) {
+            return true;
+        } else if (aX.a > bX.a) {
+            return false;
+        } else if (aX.b < bX.b) {
+            return true;
+        } else {
+            return false;
+        }
+    }
+
+    receive() external payable {}
+}
diff --git a/src/sign/BLS.sol b/src/sign/BLS.sol
index c4a8f68..66ce427 100644
--- a/src/sign/BLS.sol
+++ b/src/sign/BLS.sol
@@ -6,6 +6,10 @@ pragma solidity ^0.8.23;
 ///         defined in EIP-2537, see <https://eips.ethereum.org/EIPS/eip-2537>.
 /// @dev Precompile addresses come from the BLS addresses submodule in AlphaNet, see
 ///      <https://github.com/paradigmxyz/alphanet/blob/main/crates/precompile/src/addresses.rs>
+/// @notice `hashToCurve` logic is based on <https://github.com/ethyla/bls12-381-hash-to-curve/blob/main/src/HashToCurve.sol>
+/// with small modifications including:
+///     - Removal of low-level assembly in _modexp to ensure compatibility with EOF which does not support low-level staticcall
+///     - Usage of Fp2/G2Point structs defined here for better compatibility with existing methods
 library BLS {
     /// @dev A base field element (Fp) is encoded as 64 bytes by performing the
     /// BigEndian encoding of the corresponding (unsigned) integer. Due to the size of p,
@@ -150,4 +154,137 @@ library BLS {
         require(success, "MAP_FP2_TO_G2 failed");
         return abi.decode(output, (G2Point));
     }
+
+    /// @notice Computes a point in G2 from a message
+    /// @dev Uses the eip-2537 precompiles
+    /// @param message Arbitrarylength byte string to be hashed
+    /// @return A point in G2
+    function hashToCurveG2(bytes memory message) internal view returns (G2Point memory) {
+        // 1. u = hash_to_field(msg, 2)
+        Fp2[2] memory u = hashToFieldFp2(message, bytes("BLS_SIG_BLS12381G2_XMD:SHA-256_SSWU_RO_NUL_"));
+        // 2. Q0 = map_to_curve(u[0])
+        G2Point memory q0 = MapFp2ToG2(u[0]);
+        // 3. Q1 = map_to_curve(u[1])
+        G2Point memory q1 = MapFp2ToG2(u[1]);
+        // 4. R = Q0 + Q1
+        return G2Add(q0, q1);
+    }
+
+    /// @notice Computes a field point from a message
+    /// @dev Follows https://datatracker.ietf.org/doc/html/rfc9380#section-5.2
+    /// @param message Arbitrarylength byte string to be hashed
+    /// @param dst The domain separation tag
+    /// @return Two field points
+    function hashToFieldFp2(bytes memory message, bytes memory dst) private view returns (Fp2[2] memory) {
+        // 1. len_in_bytes = count * m * L
+        // so always 2 * 2 * 64 = 256
+        uint16 lenInBytes = 256;
+        // 2. uniform_bytes = expand_message(msg, DST, len_in_bytes)
+        bytes32[] memory pseudoRandomBytes = expandMsgXmd(message, dst, lenInBytes);
+        Fp2[2] memory u;
+        // No loop here saves 800 gas hardcoding offset an additional 300
+        // 3. for i in (0, ..., count - 1):
+        // 4.   for j in (0, ..., m - 1):
+        // 5.     elm_offset = L * (j + i * m)
+        // 6.     tv = substr(uniform_bytes, elm_offset, HTF_L)
+        // uint8 HTF_L = 64;
+        // bytes memory tv = new bytes(64);
+        // 7.     e_j = OS2IP(tv) mod p
+        // 8.   u_i = (e_0, ..., e_(m - 1))
+        // tv = bytes.concat(pseudo_random_bytes[0], pseudo_random_bytes[1]);
+        u[0].c0 = _modfield(pseudoRandomBytes[0], pseudoRandomBytes[1]);
+        u[0].c1 = _modfield(pseudoRandomBytes[2], pseudoRandomBytes[3]);
+        u[1].c0 = _modfield(pseudoRandomBytes[4], pseudoRandomBytes[5]);
+        u[1].c1 = _modfield(pseudoRandomBytes[6], pseudoRandomBytes[7]);
+        // 9. return (u_0, ..., u_(count - 1))
+        return u;
+    }
+
+    /// @notice Computes a field point from a message
+    /// @dev Follows https://datatracker.ietf.org/doc/html/rfc9380#section-5.3
+    /// @dev bytes32[] because len_in_bytes is always a multiple of 32 in our case even 128
+    /// @param message Arbitrarylength byte string to be hashed
+    /// @param dst The domain separation tag of at most 255 bytes
+    /// @param lenInBytes The length of the requested output in bytes
+    /// @return A field point
+    function expandMsgXmd(bytes memory message, bytes memory dst, uint16 lenInBytes)
+        private
+        pure
+        returns (bytes32[] memory)
+    {
+        // 1.  ell = ceil(len_in_bytes / b_in_bytes)
+        // b_in_bytes seems to be 32 for sha256
+        // ceil the division
+        uint256 ell = (lenInBytes - 1) / 32 + 1;
+
+        // 2.  ABORT if ell > 255 or len_in_bytes > 65535 or len(DST) > 255
+        require(ell <= 255, "len_in_bytes too large for sha256");
+        // Not really needed because of parameter type
+        // require(lenInBytes <= 65535, "len_in_bytes too large");
+        // no length normalizing via hashing
+        require(dst.length <= 255, "dst too long");
+
+        bytes memory dstPrime = bytes.concat(dst, bytes1(uint8(dst.length)));
+
+        // 4.  Z_pad = I2OSP(0, s_in_bytes)
+        // this should be sha256 blocksize so 64 bytes
+        bytes memory zPad = new bytes(64);
+
+        // 5.  l_i_b_str = I2OSP(len_in_bytes, 2)
+        // length in byte string?
+        bytes2 libStr = bytes2(lenInBytes);
+
+        // 6.  msg_prime = Z_pad || msg || l_i_b_str || I2OSP(0, 1) || DST_prime
+        bytes memory msgPrime = bytes.concat(zPad, message, libStr, hex"00", dstPrime);
+
+        // 7.  b_0 = H(msg_prime)
+        bytes32 b_0 = sha256(msgPrime);
+
+        bytes32[] memory b = new bytes32[](ell);
+
+        // 8.  b_1 = H(b_0 || I2OSP(1, 1) || DST_prime)
+        b[0] = sha256(bytes.concat(b_0, hex"01", dstPrime));
+
+        // 9.  for i in (2, ..., ell):
+        for (uint8 i = 2; i <= ell; i++) {
+            // 10.    b_i = H(strxor(b_0, b_(i - 1)) || I2OSP(i, 1) || DST_prime)
+            bytes memory tmp = abi.encodePacked(b_0 ^ b[i - 2], i, dstPrime);
+            b[i - 1] = sha256(tmp);
+        }
+        // 11. uniform_bytes = b_1 || ... || b_ell
+        // 12. return substr(uniform_bytes, 0, len_in_bytes)
+        // Here we don't need the uniform_bytes because b is already properly formed
+        return b;
+    }
+
+    // passing two bytes32 instead of bytes memory saves approx 700 gas per call
+    // Computes the mod against the bls12-381 field modulus
+    function _modfield(bytes32 _b1, bytes32 _b2) private view returns (Fp memory r) {
+        (bool success, bytes memory output) = address(0x5).staticcall(
+            abi.encode(
+                // arg[0] = base.length
+                0x40,
+                // arg[1] = exp.length
+                0x20,
+                // arg[2] = mod.length
+                0x40,
+                // arg[3] = base.bits
+                // places the first 32 bytes of _b1 and the last 32 bytes of _b2
+                _b1,
+                _b2,
+                // arg[4] = exp
+                // exponent always 1
+                1,
+                // arg[5] = mod
+                // this field_modulus as hex 4002409555221667393417789825735904156556882819939007885332058136124031650490837864442687629129015664037894272559787
+                // we add the 0 prefix so that the result will be exactly 64 bytes
+                // saves 300 gas per call instead of sending it along every time
+                // places the first 32 bytes and the last 32 bytes of the field modulus
+                0x000000000000000000000000000000001a0111ea397fe69a4b1ba7b6434bacd7,
+                0x64774b84f38512bf6730d2a0f6b0f6241eabfffeb153ffffb9feffffffffaaab
+            )
+        );
+        require(success, "MODEXP failed");
+        return abi.decode(output, (Fp));
+    }
 }