diff --git a/ecfactory/bls_curves/__init__.py b/ecfactory/bls_curves/__init__.py
new file mode 100644
index 0000000..aceb0db
--- /dev/null
+++ b/ecfactory/bls_curves/__init__.py
@@ -0,0 +1 @@
+from .bls_curves import make_curve
diff --git a/ecfactory/bls_curves/bls_curves.py b/ecfactory/bls_curves/bls_curves.py
new file mode 100644
index 0000000..31926e9
--- /dev/null
+++ b/ecfactory/bls_curves/bls_curves.py
@@ -0,0 +1,67 @@
+from sage.all import random, power_mod, primitive_root, Integer, random_prime, is_prime, kronecker, squarefree_part, is_square, Mod, fundamental_discriminant, sqrt, log, floor
+from ecfactory.utils import is_valid_curve
+import ecfactory.utils as utils
+import sys
+
+
+def make_curve(num_bits, num_curves=1):
+    """
+    Description:
+
+        Finds num_curves Barreto-Lynn-Scott curves with a scalar field order that is at least 2^num_bits.
+
+    Input:
+
+        num_bits - number of bits for the prime subgroup order of the curve
+        num_curves - number of curves to find
+
+    Output:
+
+        curves - list of the first num_curves BLS curves each of prime order at least 2^num_bits;
+                 each curve is represented as a tuple (q,t,r,k,D)
+
+    """
+    def R(y):
+        x = Integer(y)
+        return pow(x, 4) - pow(x, 2) + 1
+
+    def P(y):
+        x = Integer(y)
+        r = R(y)
+        return Integer(floor((pow(x-1, 2) * r) / 3 + x))
+
+    x = Integer(floor(pow(2, (num_bits)/4.0)))
+    q = 0
+    r = 0
+    t = 0
+    curve_num = 0
+    curves = []
+    while curve_num < num_curves or (log(r).n()/log(2).n() < 2*num_bits and not (utils.is_suitable_q(q) and utils.is_suitable_r(r) and utils.is_suitable_curve(q, t, r, 12, -3, num_bits))):
+        r = R(x)
+        q = P(x)
+        t = Integer(x + 1)
+
+        b = utils.is_suitable_q(q) and utils.is_suitable_r(
+            r) and utils.is_suitable_curve(q, t, r, 12, -3, num_bits)
+        if b:
+            sys.stdout.write(".")
+            try:
+                assert floor(log(r)/log(2)) + \
+                    1 >= num_bits, 'Subgroup not large enough'
+                curves.append((q, t, r, 12, -3))
+                curve_num += 1
+            except AssertionError as e:
+                pass
+        if curve_num < num_curves or not b:
+            q = P(-x)
+            t = Integer(-x + 1)
+            if (utils.is_suitable_q(q) and utils.is_suitable_r(r) and utils.is_suitable_curve(q, t, r, 12, -3, num_bits)):
+                try:
+                    assert floor(log(r)/log(2)) + \
+                        1 >= num_bits, 'Subgroup not large enough'
+                    curves.append((q, t, r, 12, -3))
+                    curve_num += 1
+                except AssertionError as e:
+                    pass
+        x += 1
+    return curves
diff --git a/ecfactory/bls_curves/bls_curves_cpp/README.md b/ecfactory/bls_curves/bls_curves_cpp/README.md
new file mode 100644
index 0000000..6f07364
--- /dev/null
+++ b/ecfactory/bls_curves/bls_curves_cpp/README.md
@@ -0,0 +1,67 @@
+# bls_curves_cpp
+
+This module provides a fast C++ implementation for finding Barreto-Lynn-Scott curve parameters.
+
+## Installation
+
+### Ubuntu 14.04
+
+Install GMP and NTL:
+
+```
+sudo apt-get install libgmp3-dev liblstl27
+```
+
+### macOS
+
+Install Homebrew:
+
+```
+ruby -e "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)" < /dev/null 2> /dev/null
+```
+
+Install GMP:
+
+```
+brew install gmp
+```
+
+Install NTL:
+
+```
+brew install ntl
+```
+
+## Run
+
+Compile with:
+
+```
+g++ -o bls_search bls_search.cpp -lntl -lgmp
+```
+
+Run with:
+
+```
+./bls_search {seed for RNG} {even_shift_of_x} {2-adicity}
+```
+
+For example, using the seed 100, and sampling values of x that are multiples of 2^2, looking for a two adicity of 10, would be:
+
+```
+./bls_search 100 2 10
+```
+
+## Log and analyze
+
+We provide a statistics script to generate a histogram of findings. To use it, log the output of `bls_search` as follows:
+
+```
+./bls_search 100 2 10 > log.txt
+```
+
+Then use `bls_statistics` on the log as follows:
+
+```
+./bls_statistics log.txt
+```
diff --git a/ecfactory/bls_curves/bls_curves_cpp/bls_search.cpp b/ecfactory/bls_curves/bls_curves_cpp/bls_search.cpp
new file mode 100644
index 0000000..431f7a8
--- /dev/null
+++ b/ecfactory/bls_curves/bls_curves_cpp/bls_search.cpp
@@ -0,0 +1,198 @@
+#include <cstdio>
+#include <iostream>
+#include <vector>
+#include <NTL/ZZ.h>
+#include <NTL/tools.h>
+
+using namespace std;
+using namespace NTL;
+
+#define REPORT_INTERVAL 100000
+#define MAX_R_BITS 255
+
+/*
+  This code modifies the method of [BLS02] to search for candidate parameters of
+  Barreto-Lynn-Scott curves with high 2-adicity. (The curve itself can be explicitly
+  constructed by following the method described in [BLS02].)
+
+  [BLS02]   = "Constructing Elliptic Curves with Prescribed Embedding Degrees"
+*/
+
+ZZ r_(ZZ x)
+{
+    ZZ x2 = x * x;
+    ZZ x3 = x2 * x;
+    ZZ x4 = x3 * x;
+    return x4 - x2 + 1;
+}
+
+ZZ q_(ZZ x)
+{
+    ZZ r = r_(x);
+    ZZ q = (x - 1) * (x - 1) * r;
+
+    bool q_is_int = q % 3 == 0;
+    if (!q_is_int)
+    {
+        return ZZ(0);
+    }
+
+    return q / 3 + x;
+}
+
+int main(int argc, char **argv)
+{
+    if (argc != 4)
+    {
+        cout << "usage: " << argv[0] << " [rand_seed even_offset wanted_two_adicity]\n";
+        return 0;
+    }
+
+    /* Collect inputs */
+    ZZ seed;
+    conv(seed, atoi(argv[1]));
+    long even_offset = atoi(argv[2]);
+    long wanted_two_adicity = atoi(argv[3]);
+
+    /* |x| ~ 64 bits so that |r| = 4 * |x| ~ 256 bits */
+    long num_x_bits = 64;
+
+    long num_iters = 0;
+    long num_found = 0;
+
+    SetSeed(seed);
+
+    /**
+     * Compute candidate BLS parameters using the formulas:
+     * r = x^4 - x^2 + 1
+     * q = ((x-1)^2 * r) / 3 + x
+     * (see [BLS02])
+     */
+
+    while (1)
+    {
+        num_iters += 4;
+
+        /* Sample x */
+        ZZ x = RandomLen_ZZ(num_x_bits);
+
+        /**
+         * Make x even and divisible by a large power of 2 to improve two adicity.
+         * The resulting r is s.t. -1 is a square in Fq.
+         */
+
+        x = ((x >> even_offset) << even_offset);
+        ZZ r = r_(x);
+        ZZ q = q_(x);
+
+        long num_q_bits = NumBits(q);
+        long num_r_bits = NumBits(r);
+        long two_adicity = NumTwos(r - 1);
+
+        if (num_r_bits > MAX_R_BITS)
+        {
+            continue;
+        }
+
+        if (ProbPrime(r) && ProbPrime(q) && (two_adicity >= wanted_two_adicity))
+        {
+            cout << "x = " << x << "\n";
+            cout << "q = " << q << "\n";
+            cout << "r = " << r << "\n";
+            cout << "log2(q) =" << num_q_bits << "\n";
+            cout << "log2(r) =" << num_r_bits << "\n";
+            cout << "ord_2(r-1) = " << two_adicity << "\n";
+            cout.flush();
+            ++num_found;
+        }
+
+        /**
+         * -x will yield the same r, but a different q.
+         */
+
+        q = q_(-x);
+
+        num_q_bits = NumBits(q);
+        num_r_bits = NumBits(r);
+        two_adicity = NumTwos(r - 1);
+
+        if (num_r_bits > MAX_R_BITS)
+        {
+            continue;
+        }
+
+        if (ProbPrime(r) && ProbPrime(q) && (two_adicity >= wanted_two_adicity))
+        {
+            cout << "x = " << -x << "\n";
+            cout << "q = " << q << "\n";
+            cout << "r = " << r << "\n";
+            cout << "log2(q) =" << num_q_bits << "\n";
+            cout << "log2(r) =" << num_r_bits << "\n";
+            cout << "ord_2(r-1) = " << two_adicity << "\n";
+            cout.flush();
+            ++num_found;
+        }
+
+        /**
+         * An x of the form k.2^n + 1 also yields a scalar field with high (n+1) two adicity.
+         */
+
+        x = x + 1;
+        r = r_(x);
+        q = q_(x);
+
+        num_q_bits = NumBits(q);
+        num_r_bits = NumBits(r);
+        two_adicity = NumTwos(r - 1);
+
+        if (num_r_bits > MAX_R_BITS)
+        {
+            continue;
+        }
+
+        if (ProbPrime(r) && ProbPrime(q) && (two_adicity >= wanted_two_adicity))
+        {
+            cout << "x = " << x << "\n";
+            cout << "q = " << q << "\n";
+            cout << "r = " << r << "\n";
+            cout << "log2(q) =" << num_q_bits << "\n";
+            cout << "log2(r) =" << num_r_bits << "\n";
+            cout << "ord_2(r-1) = " << two_adicity << "\n";
+            cout.flush();
+            ++num_found;
+        }
+
+        /**
+         * -x will yield the same r, but a different q.
+         */
+
+        q = q_(-x);
+
+        num_q_bits = NumBits(q);
+        num_r_bits = NumBits(r);
+        two_adicity = NumTwos(r - 1);
+
+        if (num_r_bits > MAX_R_BITS)
+        {
+            continue;
+        }
+
+        if (ProbPrime(r) && ProbPrime(q) && (two_adicity >= wanted_two_adicity))
+        {
+            cout << "x = " << -x << "\n";
+            cout << "q = " << q << "\n";
+            cout << "r = " << r << "\n";
+            cout << "log2(q) =" << num_q_bits << "\n";
+            cout << "log2(r) =" << num_r_bits << "\n";
+            cout << "ord_2(r-1) = " << two_adicity << "\n";
+            cout.flush();
+            ++num_found;
+        }
+
+        if (num_iters % REPORT_INTERVAL == 0)
+        {
+            printf("[ num_iters = %ld , num_found = %0.2f per %d ]\n", num_iters, 1. * REPORT_INTERVAL * num_found / num_iters, REPORT_INTERVAL);
+            fflush(stdout);
+        }
+    }
+}
diff --git a/ecfactory/bls_curves/bls_curves_cpp/bls_statistics b/ecfactory/bls_curves/bls_curves_cpp/bls_statistics
new file mode 100755
index 0000000..170b910
--- /dev/null
+++ b/ecfactory/bls_curves/bls_curves_cpp/bls_statistics
@@ -0,0 +1,9 @@
+#!/bin/bash
+
+FILE=$1
+
+for i in `seq 10 64`;
+do
+	echo "$i:"
+	cat $FILE | grep "(r-1) = $i" | wc -l
+done
diff --git a/ecfactory/bls_curves/bls_curves_examples.py b/ecfactory/bls_curves/bls_curves_examples.py
new file mode 100644
index 0000000..2edaa28
--- /dev/null
+++ b/ecfactory/bls_curves/bls_curves_examples.py
@@ -0,0 +1,9 @@
+import ecfactory.bls_curves as bls
+from ecfactory.utils import print_curve
+
+# Example
+num_bits = 100  # number of bits in r
+num_curves = 10  # number of curves to find
+curves = bls.make_curve(num_bits, num_curves)
+for q, t, r, k, D in curves:
+    print_curve(q, t, r, k, D)
diff --git a/ecfactory/bls_curves/bls_curves_tests.py b/ecfactory/bls_curves/bls_curves_tests.py
new file mode 100644
index 0000000..d6ef164
--- /dev/null
+++ b/ecfactory/bls_curves/bls_curves_tests.py
@@ -0,0 +1,29 @@
+from sage.all import randint
+import ecfactory.bls_curves as bls
+import ecfactory.utils as utils
+
+
+def test_bls(num_tests):
+    # finds num_tests BLS curves with random bit sizes;
+    # checks that each curve found is a suitable curve
+    print('testing BLS...')
+    fail = 0
+    for i in range(0, num_tests):
+        num_bits = randint(50, 100)
+        num_curves = randint(1, 20)
+        try:
+            curves = bls.make_curve(num_bits, num_curves)
+            assert len(curves) == num_curves
+            for q, t, r, k, D in curves:
+                assert utils.is_suitable_curve(q, t, r, k, D, num_bits)
+        except AssertionError as e:
+            fail += 1
+    if fail == 0:
+        print('test passed')
+        return True
+    else:
+        print("failed %.2f" % (100*fail/num_tests) + "% of tests!")
+        return False
+
+
+test_bls(10)
diff --git a/ecfactory/bls_curves/readme.md b/ecfactory/bls_curves/readme.md
new file mode 100644
index 0000000..86c1ef3
--- /dev/null
+++ b/ecfactory/bls_curves/readme.md
@@ -0,0 +1,41 @@
+Barreto-Lynn-Scott Curves
+======================
+
+Overview
+--------
+
+This module provides functionality to construct _Barreto-Lynn-Scott (BLS) curves_, following the procedure described in [\[BLS02\]](/references/Barreto%20Lynn%20Scott%202002%20---%20Constructing%20Elliptic%20Curves%20with%20Prescribed%20Embedding%20Degrees.pdf). The module contains three files:
+
+* `bls_curves.py`, which contains the algorithm to construct Barreto-Lynn-Scott curves;
+
+* `bls_curves_examples.py`, which contains code examples;
+
+* `bls_curves_tests.py`, which contains unit tests.
+
+Throughout,
+_q_ denotes the prime size of the base field;
+_t_ denotes the trace of Frobenius;
+_r_ denotes the prime size of the group;
+_k_ denotes the embedding degree;
+_D_ denotes the (negative) fundamental discriminant.
+
+Methods
+-------
+
+The main method in `bls_curves.py` is
+
+```python
+make_curve(num_bits, num_curves=1)
+```
+
+It outputs a list of the first *num\_curves* BLS curves of prime order at least 2<sup>*num\_bits*</sup>.
+
+Examples
+--------
+
+The code example in `bls_curves_examples.py` shows how to find 10 BLS curves with a prime order that is at least 2<sup>100</sup>.
+
+Tests
+-----
+
+The test in `bls_curves_tests.py` runs the algorithm on random bit-sizes and checks validity of the output.
diff --git a/references/Barreto Lynn Scott 2002 --- Constructing Elliptic Curves with Prescribed Embedding Degrees.pdf b/references/Barreto Lynn Scott 2002 --- Constructing Elliptic Curves with Prescribed Embedding Degrees.pdf
new file mode 100644
index 0000000..687e3fa
Binary files /dev/null and b/references/Barreto Lynn Scott 2002 --- Constructing Elliptic Curves with Prescribed Embedding Degrees.pdf differ