diff --git a/eth/constants.py b/eth/constants.py
index 4388e35ac9..4a826b7ddd 100644
--- a/eth/constants.py
+++ b/eth/constants.py
@@ -100,6 +100,15 @@
 GAS_ECPAIRING_BASE = 100000
 GAS_ECPAIRING_PER_POINT = 80000
 
+GAS_BLS_G1_ADD = 600
+GAS_BLS_G1_MUL = 12000
+GAS_BLS_G2_ADD = 4500
+GAS_BLS_G2_MUL = 55000
+GAS_BLS_PAIRING_BASE = 115000
+GAS_BLS_PAIRING_PER_PAIR = 23000
+GAS_BLS_MAP_FP_TO_G1 = 5500
+GAS_BLS_MAP_FP2_TO_G2 = 110000
+
 
 #
 # Gas Limit
diff --git a/eth/precompiles/__init__.py b/eth/precompiles/__init__.py
index e18c7e95f7..7a7b2d34d4 100644
--- a/eth/precompiles/__init__.py
+++ b/eth/precompiles/__init__.py
@@ -7,3 +7,14 @@
 from .ecmul import ecmul  # noqa: F401
 from .ecpairing import ecpairing  # noqa: F401
 from .blake2 import blake2b_fcompress  # noqa: F401
+from .bls import (  # noqa: F401
+    g1_add as bls_g1_add,
+    g1_mul as bls_g1_mul,
+    g1_multiexp as bls_g1_multiexp,
+    g2_add as bls_g2_add,
+    g2_mul as bls_g2_mul,
+    g2_multiexp as bls_g2_multiexp,
+    pairing as bls_pairing,
+    map_fp_to_g1 as bls_map_fp_to_g1,
+    map_fp2_to_g2 as bls_map_fp2_to_g2,
+)
diff --git a/eth/precompiles/bls.py b/eth/precompiles/bls.py
new file mode 100644
index 0000000000..fdc875b06b
--- /dev/null
+++ b/eth/precompiles/bls.py
@@ -0,0 +1,236 @@
+from typing import Tuple
+
+from eth_utils import (
+    ValidationError,
+    big_endian_to_int,
+)
+from py_ecc import (
+    optimized_bls12_381 as bls12_381,
+    bls
+)
+
+from eth import constants
+from eth.exceptions import (
+    VMError,
+)
+
+from eth.vm.computation import (
+    BaseComputation,
+)
+
+
+FP2_SIZE_IN_BYTES = 128
+G1_SIZE_IN_BYTES = 128
+G2_SIZE_IN_BYTES = 256
+# Constant for parsing pairs of points during pairing check
+G1_TO_G2_OFFSET = G1_SIZE_IN_BYTES + G2_SIZE_IN_BYTES
+
+G1Point = Tuple[bls12_381.FQ, bls12_381.FQ]
+G2Point = Tuple[bls12_381.FQ2, bls12_381.FQ2]
+
+
+def _parse_g1_point(data: bytes) -> G1Point:
+    if len(data) != G1_SIZE_IN_BYTES:
+        raise ValidationError("invalid size of G1 input")
+
+    x = bls12_381.FQ(int.from_bytes(data[0:64], byteorder="big"))
+    y = bls12_381.FQ(int.from_bytes(data[64:128], byteorder="big"))
+    point = (x, y)
+
+    if not bls12_381.is_on_curve((x, y, bls12_381.FQ.one()), bls12_381.b):
+        raise ValidationError("invalid G1 point not on curve")
+
+    return point
+
+
+def _serialize_g1(result: G1Point) -> bytes:
+    return b"".join(
+        (
+            int(result[0]).to_bytes(64, byteorder="big"),
+            int(result[1]).to_bytes(64, byteorder="big"),
+        )
+    )
+
+
+def g1_add(computation: BaseComputation,
+           gas_cost: int = constants.GAS_BLS_G1_ADD) -> BaseComputation:
+    raise NotImplementedError()
+
+
+def g1_mul(computation: BaseComputation,
+           gas_cost: int = constants.GAS_BLS_G1_MUL) -> BaseComputation:
+    raise NotImplementedError()
+
+
+def g1_multiexp(computation: BaseComputation) -> BaseComputation:
+    # NOTE: gas cost involves a discount based on the number of points involved
+    # TODO load discount table and compute gas cost based on number of inputs
+    raise NotImplementedError()
+
+
+def _parse_g2_point(data: bytes) -> G2Point:
+    if len(data) != G2_SIZE_IN_BYTES:
+        raise ValidationError("invalid size of G2 input")
+
+    x = bls12_381.FQ2(
+        (
+            int.from_bytes(data[0:64], byteorder="big"),
+            int.from_bytes(data[64:128], byteorder="big")
+        )
+    )
+    y = bls12_381.FQ2(
+        (
+            int.from_bytes(data[128:192], byteorder="big"),
+            int.from_bytes(data[192:256], byteorder="big")
+        )
+    )
+    point = (x, y)
+
+    if not bls12_381.is_on_curve((x, y, bls12_381.FQ2.one()), bls12_381.b2):
+        raise ValidationError("invalid G2 point not on curve")
+
+    return point
+
+
+def _serialize_g2(result: G2Point) -> bytes:
+    return b"".join(
+        (
+            result[0].coeffs[0].to_bytes(64, byteorder="big"),
+            result[0].coeffs[1].to_bytes(64, byteorder="big"),
+            result[1].coeffs[0].to_bytes(64, byteorder="big"),
+            result[1].coeffs[1].to_bytes(64, byteorder="big"),
+        )
+    )
+
+
+def _g2_add(x: G2Point, y: G2Point) -> G2Point:
+    result = bls12_381.add((x[0], x[1], bls12_381.FQ2.one()), (y[0], y[1], bls12_381.FQ2.one()))
+    return bls12_381.normalize(result)
+
+
+def g2_add(computation: BaseComputation,
+           gas_cost: int = constants.GAS_BLS_G2_ADD) -> BaseComputation:
+    computation.consume_gas(gas_cost, reason='BLS_G2_ADD Precompile')
+
+    try:
+        input_data = computation.msg.data_as_bytes
+        x = _parse_g2_point(input_data[:G2_SIZE_IN_BYTES])
+        y = _parse_g2_point(input_data[G2_SIZE_IN_BYTES:])
+        result = _g2_add(x, y)
+    except ValidationError:
+        raise VMError("Invalid BLS_G2_ADD parameters")
+
+    computation.output = _serialize_g2(result)
+    return computation
+
+
+def _g2_mul(x: G2Point, k: int) -> G2Point:
+    result = bls12_381.multiply((x[0], x[1], bls12_381.FQ2.one()), k)
+    return bls12_381.normalize(result)
+
+
+def _parse_scalar(data: bytes) -> int:
+    if len(data) != 32:
+        raise ValidationError("invalid size of scalar input")
+
+    return big_endian_to_int(data)
+
+
+def g2_mul(computation: BaseComputation,
+           gas_cost: int = constants.GAS_BLS_G2_MUL) -> BaseComputation:
+    computation.consume_gas(gas_cost, reason='BLS_G2_MUL Precompile')
+
+    try:
+        input_data = computation.msg.data_as_bytes
+        x = _parse_g2_point(input_data[:G2_SIZE_IN_BYTES])
+        k = _parse_scalar(input_data[G2_SIZE_IN_BYTES:])
+        result = _g2_mul(x, k)
+    except ValidationError:
+        raise VMError("Invalid BLS_G2_MUL parameters")
+
+    computation.output = _serialize_g2(result)
+    return computation
+
+
+def g2_multiexp(computation: BaseComputation) -> BaseComputation:
+    # NOTE: gas cost involves a discount based on the number of points involved
+    # TODO load discount table and compute gas cost based on number of inputs
+    raise NotImplementedError()
+
+
+def _pairing(input_data: bytes) -> bool:
+    field_element = bls12_381.FQ12.one()
+    for next_index in range(0, len(input_data), 384):
+        p = _parse_g1_point(input_data[next_index:next_index + G1_SIZE_IN_BYTES])
+        q = _parse_g2_point(
+            input_data[next_index + G1_SIZE_IN_BYTES:next_index + G1_TO_G2_OFFSET]
+        )
+        projective_p = (p[0], p[1], bls12_381.FQ.one())
+        projective_q = (q[0], q[1], bls12_381.FQ2.one())
+        field_element *= bls12_381.pairing(projective_q, projective_p, final_exponentiate=False)
+
+    return bls12_381.final_exponentiate(field_element) == bls12_381.FQ12.one()
+
+
+def _serialize_boolean(value: bool) -> bytes:
+    return int(value).to_bytes(32, byteorder="big")
+
+
+def pairing(computation: BaseComputation,
+            gas_cost_base: int = constants.GAS_BLS_PAIRING_BASE,
+            gas_cost_per_pair: int = constants.GAS_BLS_PAIRING_PER_PAIR) -> BaseComputation:
+    input_data = computation.msg.data_as_bytes
+    if len(input_data) % 384:
+        # data length must be an exact multiple of 384
+        raise VMError("Invalid BLS_PAIRING parameters")
+
+    num_points = len(input_data) // 384
+    gas_cost = gas_cost_base + num_points * gas_cost_per_pair
+
+    computation.consume_gas(gas_cost, reason='BLS_PAIRING Precompile')
+
+    try:
+        result = _pairing(input_data)
+    except ValidationError:
+        raise VMError("Invalid BLS_PAIRING parameters")
+
+    computation.output = _serialize_boolean(result)
+    return computation
+
+
+def map_fp_to_g1(computation: BaseComputation,
+                 gas_cost: int = constants.GAS_BLS_MAP_FP_TO_G1) -> BaseComputation:
+    raise NotImplementedError()
+
+
+def _parse_fp2_element(data: bytes) -> bls12_381.FQ2:
+    if len(data) != FP2_SIZE_IN_BYTES:
+        raise ValidationError("invalid size of FP2 input")
+
+    return bls12_381.FQ2(
+        (
+            int.from_bytes(data[:64], byteorder="big"),
+            int.from_bytes(data[64:], byteorder="big")
+        )
+    )
+
+
+def _map_fp2_to_g2(field_element: bls12_381.FQ2) -> G2Point:
+    point = bls.hash_to_curve.map_to_curve_G2(field_element)
+    group_element = bls.hash_to_curve.clear_cofactor_G2(point)
+    return bls12_381.normalize(group_element)
+
+
+def map_fp2_to_g2(computation: BaseComputation,
+                  gas_cost: int = constants.GAS_BLS_MAP_FP2_TO_G2) -> BaseComputation:
+    computation.consume_gas(gas_cost, reason='BLS_MAP_FP2_TO_G2 Precompile')
+
+    try:
+        input_data = computation.msg.data_as_bytes
+        field_element = _parse_fp2_element(input_data[:FP2_SIZE_IN_BYTES])
+        result = _map_fp2_to_g2(field_element)
+    except ValidationError:
+        raise VMError("Invalid BLS_MAP_FP2_TO_G2 parameters")
+
+    computation.output = _serialize_g2(result)
+    return computation
diff --git a/eth/vm/forks/berlin/__init__.py b/eth/vm/forks/berlin/__init__.py
new file mode 100644
index 0000000000..4e521f76c0
--- /dev/null
+++ b/eth/vm/forks/berlin/__init__.py
@@ -0,0 +1,31 @@
+from typing import (
+    Type,
+)
+
+from eth.rlp.blocks import BaseBlock
+from eth.vm.forks.constantinople import (
+    ConstantinopleVM,
+)
+from eth.vm.state import BaseState
+
+from .blocks import BerlinBlock
+from .headers import (
+    compute_berlin_difficulty,
+    configure_berlin_header,
+    create_berlin_header_from_parent,
+)
+from .state import BerlinState
+
+
+class BerlinVM(ConstantinopleVM):
+    # fork name
+    fork = 'berlin'
+
+    # classes
+    block_class: Type[BaseBlock] = BerlinBlock
+    _state_class: Type[BaseState] = BerlinState
+
+    # Methods
+    create_header_from_parent = staticmethod(create_berlin_header_from_parent)  # type: ignore
+    compute_difficulty = staticmethod(compute_berlin_difficulty)    # type: ignore
+    configure_header = configure_berlin_header
diff --git a/eth/vm/forks/berlin/blocks.py b/eth/vm/forks/berlin/blocks.py
new file mode 100644
index 0000000000..61df903a39
--- /dev/null
+++ b/eth/vm/forks/berlin/blocks.py
@@ -0,0 +1,22 @@
+from rlp.sedes import (
+    CountableList,
+)
+from eth.rlp.headers import (
+    BlockHeader,
+)
+from eth.vm.forks.muir_glacier.blocks import (
+    MuirGlacierBlock,
+)
+
+from .transactions import (
+    BerlinTransaction,
+)
+
+
+class BerlinBlock(MuirGlacierBlock):
+    transaction_class = BerlinTransaction
+    fields = [
+        ('header', BlockHeader),
+        ('transactions', CountableList(transaction_class)),
+        ('uncles', CountableList(BlockHeader))
+    ]
diff --git a/eth/vm/forks/berlin/computation.py b/eth/vm/forks/berlin/computation.py
new file mode 100644
index 0000000000..ecffde5919
--- /dev/null
+++ b/eth/vm/forks/berlin/computation.py
@@ -0,0 +1,41 @@
+from eth_utils.toolz import (
+    merge,
+)
+
+from eth import precompiles
+from eth._utils.address import (
+    force_bytes_to_address,
+)
+from eth.vm.forks.muir_glacier.computation import (
+    MUIR_GLACIER_PRECOMPILES
+)
+from eth.vm.forks.muir_glacier.computation import (
+    MuirGlacierComputation,
+)
+
+from .opcodes import BERLIN_OPCODES
+
+BERLIN_PRECOMPILES = merge(
+    MUIR_GLACIER_PRECOMPILES,
+    {
+        force_bytes_to_address(b'\x0a'): precompiles.bls_g1_add,
+        force_bytes_to_address(b'\x0b'): precompiles.bls_g1_mul,
+        force_bytes_to_address(b'\x0c'): precompiles.bls_g1_multiexp,
+        force_bytes_to_address(b'\x0d'): precompiles.bls_g2_add,
+        force_bytes_to_address(b'\x0e'): precompiles.bls_g2_mul,
+        force_bytes_to_address(b'\x0f'): precompiles.bls_g2_multiexp,
+        force_bytes_to_address(b'\x10'): precompiles.bls_pairing,
+        force_bytes_to_address(b'\x11'): precompiles.bls_map_fp_to_g1,
+        force_bytes_to_address(b'\x12'): precompiles.bls_map_fp2_to_g2,
+    }
+)
+
+
+class BerlinComputation(MuirGlacierComputation):
+    """
+    A class for all execution computations in the ``Berlin`` fork.
+    Inherits from :class:`~eth.vm.forks.muir_glacier.MuirGlacierComputation`
+    """
+    # Override
+    opcodes = BERLIN_OPCODES
+    _precompiles = BERLIN_PRECOMPILES
diff --git a/eth/vm/forks/berlin/headers.py b/eth/vm/forks/berlin/headers.py
new file mode 100644
index 0000000000..09352021f2
--- /dev/null
+++ b/eth/vm/forks/berlin/headers.py
@@ -0,0 +1,13 @@
+from eth.vm.forks.muir_glacier.headers import (
+    configure_header,
+    create_header_from_parent,
+    compute_muir_glacier_difficulty,
+)
+
+
+compute_berlin_difficulty = compute_muir_glacier_difficulty
+
+create_berlin_header_from_parent = create_header_from_parent(
+    compute_berlin_difficulty
+)
+configure_berlin_header = configure_header(compute_berlin_difficulty)
diff --git a/eth/vm/forks/berlin/opcodes.py b/eth/vm/forks/berlin/opcodes.py
new file mode 100644
index 0000000000..a0b7e2bd21
--- /dev/null
+++ b/eth/vm/forks/berlin/opcodes.py
@@ -0,0 +1,18 @@
+import copy
+
+from eth_utils.toolz import merge
+
+
+from eth.vm.forks.muir_glacier.opcodes import (
+    MUIR_GLACIER_OPCODES,
+)
+
+
+UPDATED_OPCODES = {
+    # New opcodes
+}
+
+BERLIN_OPCODES = merge(
+    copy.deepcopy(MUIR_GLACIER_OPCODES),
+    UPDATED_OPCODES,
+)
diff --git a/eth/vm/forks/berlin/state.py b/eth/vm/forks/berlin/state.py
new file mode 100644
index 0000000000..c5e9d85ed5
--- /dev/null
+++ b/eth/vm/forks/berlin/state.py
@@ -0,0 +1,9 @@
+from eth.vm.forks.muir_glacier.state import (
+    MuirGlacierState
+)
+
+from .computation import BerlinComputation
+
+
+class BerlinState(MuirGlacierState):
+    computation_class = BerlinComputation
diff --git a/eth/vm/forks/berlin/transactions.py b/eth/vm/forks/berlin/transactions.py
new file mode 100644
index 0000000000..548a794a7c
--- /dev/null
+++ b/eth/vm/forks/berlin/transactions.py
@@ -0,0 +1,42 @@
+from eth_keys.datatypes import PrivateKey
+from eth_typing import Address
+
+from eth.vm.forks.muir_glacier.transactions import (
+    MuirGlacierTransaction,
+    MuirGlacierUnsignedTransaction,
+)
+
+from eth._utils.transactions import (
+    create_transaction_signature,
+)
+
+
+class BerlinTransaction(MuirGlacierTransaction):
+    @classmethod
+    def create_unsigned_transaction(cls,
+                                    *,
+                                    nonce: int,
+                                    gas_price: int,
+                                    gas: int,
+                                    to: Address,
+                                    value: int,
+                                    data: bytes) -> 'BerlinUnsignedTransaction':
+        return BerlinUnsignedTransaction(nonce, gas_price, gas, to, value, data)
+
+
+class BerlinUnsignedTransaction(MuirGlacierUnsignedTransaction):
+    def as_signed_transaction(self,
+                              private_key: PrivateKey,
+                              chain_id: int=None) -> BerlinTransaction:
+        v, r, s = create_transaction_signature(self, private_key, chain_id=chain_id)
+        return BerlinTransaction(
+            nonce=self.nonce,
+            gas_price=self.gas_price,
+            gas=self.gas,
+            to=self.to,
+            value=self.value,
+            data=self.data,
+            v=v,
+            r=r,
+            s=s,
+        )
diff --git a/tests/core/precompiles/test_bls.py b/tests/core/precompiles/test_bls.py
new file mode 100644
index 0000000000..0df5bad1df
--- /dev/null
+++ b/tests/core/precompiles/test_bls.py
@@ -0,0 +1,10 @@
+from py_ecc.optimized_bls12_381.optimized_curve import G1, G2, neg, FQ12
+from eth.precompiles.bls import _pairing, _serialize_g2, _serialize_g1
+
+def test_pairing_precompile():
+    # assert pairing(G1, G2) * pairing(neg(G1), G2) == FQ12.one()
+    serialized_G1 = _serialize_g1(G1)
+    serialized_G2 = _serialize_g2(G2)
+    serialized_neg_G1 = _serialize_g1(neg(G1))
+    input_data = serialized_G1 + serialized_G2 + serialized_neg_G1 + serialized_G2
+    assert _pairing(input_data)