diff --git a/LICENSE b/LICENSE
new file mode 100644
index 0000000..59312e9
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2021 Nicolas Tessore
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
diff --git a/MANIFEST.in b/MANIFEST.in
new file mode 100644
index 0000000..e7f8c58
--- /dev/null
+++ b/MANIFEST.in
@@ -0,0 +1,3 @@
+include flt.pyx
+include dctdlt.c
+exclude flt.c
diff --git a/README.md b/README.md
new file mode 100644
index 0000000..fc4191b
--- /dev/null
+++ b/README.md
@@ -0,0 +1,20 @@
+
+flt
+===
+
+**fast Legendre transform**
+
+This is a minimal Python package for fast discrete Legendre transforms (DLTs).
+The implementation uses a recursive version of the matrix relations by Alpert &
+Rokhlin (1991) to compute the DLT via a discrete cosine transform (DCT).
+
+The package can be installed using pip:
+
+    pip install flt
+
+For more information, please see the [documentation].
+
+Current functionality covers the absolutely minimal use case.  Please open an
+issue on GitHub if you would like to see anything added.
+
+[documentation]: https://cltools.readthedocs.io/flt/
diff --git a/dctdlt.c b/dctdlt.c
new file mode 100644
index 0000000..2ab2901
--- /dev/null
+++ b/dctdlt.c
@@ -0,0 +1,107 @@
+// dctdlt.c
+// ========
+// discrete Legendre transform via DCT
+//
+// author: Nicolas Tessore <n.tessore@ucl.ac.uk>
+// license: MIT
+//
+// Synopsis
+// --------
+// The `dctdlt` and `dltdct` functions convert the coefficients of a discrete
+// cosine transform (DCT) to the coefficients of a discrete Legendre transform
+// (DLT) and vice versa [1].
+//
+// References
+// ----------
+// [1] Alpert, B. K., & Rokhlin, V. (1991). A fast algorithm for the evaluation
+//     of Legendre expansions. SIAM Journal on Scientific and Statistical
+//     Computing, 12(1), 158-179.
+//
+
+#define DCTDLT_VERSION 20210318L
+
+
+// dctdlt
+// ======
+// convert DCT coefficients to DLT coefficients
+//
+// Parameters
+// ----------
+// n : unsigned int
+//     Length of the input array.
+// dct : (n,) array of double
+//     Input DCT coefficients.
+// dlt : (n,) array of double, output
+//     Output DLT coefficients.
+//
+void dctdlt(unsigned int n, const double* dct, double* dlt)
+{
+    double a, b;
+    unsigned int k, l;
+
+    // first row
+    a = 1.;
+    b = a;
+    dlt[0] = 0.5*b*dct[0];
+    for(k = 2; k < n; k += 2)
+    {
+        b *= (k-3.)/(k+1.);
+        dlt[0] += b*dct[k];
+    }
+
+    // remaining rows
+    for(l = 1; l < n; ++l)
+    {
+        a /= (1. - 0.5/l);
+        b = a;
+        dlt[l] = b*dct[l];
+        for(k = l+2; k < n; k += 2)
+        {
+            b *= (k*(k+l-2.)*(k-l-3.))/((k-2.)*(k+l+1.)*(k-l));
+            dlt[l] += b*dct[k];
+        }
+    }
+}
+
+
+// dltdct
+// ======
+// convert DLT coefficients to DCT coefficients
+//
+// Parameters
+// ----------
+// n : unsigned int
+//     Length of the input array.
+// dlt : (n,) array of double
+//     Input DLT coefficients.
+// dct : (n,) array of double, output
+//     Output DCT coefficients.
+//
+void dltdct(unsigned int n, const double* dlt, double* dct)
+{
+    double a, b;
+    unsigned int k, l;
+
+    // first row
+    a = 1.;
+    b = a;
+    dct[0] = b*dlt[0];
+    for(l = 2; l < n; l += 2)
+    {
+        b *= ((l-1.)*(l-1.))/(l*l);
+        dct[0] += b*dlt[l];
+    }
+
+    // remaining rows
+    for(k = 1; k < n; ++k)
+    {
+        a *= (1. - 0.5/k);
+        b = a;
+        dct[k] = b*dlt[k];
+        for(l = k+2; l < n; l += 2)
+        {
+            b *= ((l-k-1.)*(l+k-1.))/((l-k)*(l+k));
+            dct[k] += b*dlt[l];
+        }
+    }
+}
diff --git a/flt.pyx b/flt.pyx
new file mode 100644
index 0000000..426455f
--- /dev/null
+++ b/flt.pyx
@@ -0,0 +1,377 @@
+# flt: fast Legendre transform
+#
+# author: Nicolas Tessore <n.tessore@ucl.ac.uk>
+# license: MIT
+#
+# cython: language_level=3, boundscheck=False, embedsignature=True
+#
+'''
+
+Discrete Legendre Transform (:mod:`flt`)
+========================================
+
+This is a minimal Python package for fast discrete Legendre transforms (DLTs).
+The implementation uses a recursive version of the matrix relations by Alpert &
+Rokhlin (1991) to compute the DLT via a discrete cosine transform (DCT).
+
+The package can be installed using pip::
+
+    pip install flt
+
+Current functionality covers the absolutely minimal use case.  Please open an
+issue on GitHub if you would like to see anything added.
+
+
+Reference/API
+-------------
+
+.. autosummary::
+   :toctree: api
+   :nosignatures:
+
+   dlt
+   idlt
+   dltmtx
+   idltmtx
+   theta
+
+'''
+
+__all__ = [
+    'dlt',
+    'idlt',
+    'dltmtx',
+    'idltmtx',
+    'theta',
+]
+
+
+import numpy as np
+from scipy.fft import dct, idct
+
+cdef extern from "dctdlt.c":
+    void dctdlt(unsigned int, const double*, double*)
+    void dltdct(unsigned int, const double*, double*)
+
+
+def dlt(a, closed=False):
+    r'''discrete Legendre transform
+
+    Takes a function in :math:`\mathtt{n}` points and returns its discrete
+    Legendre transform coefficients from :math:`0` to :math:`\mathtt{lmax}
+    = \mathtt{n}-1`.
+
+    The function must be given at the points :math:`\cos(\theta)` returned by
+    :func:`flt.theta`.  These can be distributed either over the open interval
+    :math:`\theta \in (0, \pi)`, or over the closed interval :math:`\theta \in
+    [0, \pi]`, in which case :math:`\theta_0 = 0`\ and :math:`\theta_{n-1} =
+    \pi`.
+
+    Parameters
+    ----------
+    a : (n,) array_like
+        Function values.
+    closed : bool, optional
+        Compute DLT over open (``closed=False``) or closed (``closed=True``)
+        interval.
+
+    Returns
+    -------
+    b : (n,) array_like
+        Legendre coefficients :math:`0` to :math:`\mathtt{lmax}`.
+
+    See Also
+    --------
+    flt.idlt : the inverse operation
+    flt.theta : compute the angles at which the function is evaluated
+
+    Notes
+    -----
+    The discrete Legendre transform takes a function :math:`f(\cos\theta)` over
+    the domain :math:`\theta \in [0, \pi]` and returns the coefficients of the
+    series
+
+    .. math::
+
+        f(\cos\theta) = \sum_{l=0}^{l_{\max}} c_l \, P_l(\cos\theta) \;,
+
+    where :math:`P_l(\cos\theta)` is a Legendre polynomial.
+
+    The computation is done in two steps:
+
+    First, the function is transformed to the coefficients of a discrete cosine
+    transform (DCT) using an inverse DCT-III for the open interval, or an
+    inverse DCT-I for the closed interval.
+
+    Second, the DCT coefficients are transformed to the DLT coefficientsu using
+    a recursive version of the matrix relation given by [1]_.
+
+    References
+    ----------
+    .. [1] Alpert, B. K., & Rokhlin, V. (1991). A fast algorithm for the
+       evaluation of Legendre expansions. SIAM Journal on Scientific and
+       Statistical Computing, 12(1), 158-179.
+
+    '''
+
+    # length n of the transform
+    if np.ndim(a) != 1:
+        raise TypeError('array must be 1d')
+    n = np.shape(a)[-1]
+
+    # type of the DCT depends on open or closed interval
+    if closed:
+        dcttype = 1
+    else:
+        dcttype = 3
+
+    # compute the DCT coefficients
+    a = idct(a, type=dcttype, axis=-1, norm=None)
+
+    # this holds the DLT coefficients
+    b = np.empty(n, dtype=float)
+
+    # these are memviews on a and b for C interop
+    cdef double[::1] a_ = a
+    cdef double[::1] b_ = b
+
+    # transform DCT coefficients to DLT coefficients using C function
+    dctdlt(n, &a_[0], &b_[0])
+
+    # done
+    return b
+
+
+def idlt(b, closed=False):
+    r'''inverse discrete Legendre transform
+
+    Takes the :math:`\mathtt{n} = \mathtt{lmax}+1` coefficients of a DLT and
+    returns the corresponding function in :math:`\mathtt{n}` points.
+
+    The function will be given at the points :math:`\cos(\theta)` returned by
+    :func:`flt.theta`.  These can be distributed either over the open interval
+    :math:`\theta \in (0, \pi)`, or over the closed interval :math:`\theta \in
+    [0, \pi]`, in which case :math:`\theta_0 = 0` and :math:`\theta_{n-1} =
+    \pi`.
+
+    Parameters
+    ----------
+    b : (n,) array_like
+        DLT coefficients from :math:`0` to :math:`\mathtt{lmax}`.
+    closed : bool, optional
+        Compute function over open (``closed=False``) or closed
+        (``closed=True``) interval.
+
+    Returns
+    -------
+    a : (n,) array_like
+        Function values.
+
+    See Also
+    --------
+    flt.dlt : the forward operation
+    flt.theta : compute the angles at which the function is evaluated
+
+    Notes
+    -----
+    The inverse discrete Legendre transform returns a function
+    :math:`f(\cos\theta)` over the domain :math:`\theta \in [0, \pi]` given the
+    coefficients of the series
+
+    .. math::
+
+        f(\cos\theta) = \sum_{l=0}^{l_{\max}} c_l \, P_l(\cos\theta) \;,
+
+    where :math:`P_l(\cos\theta)` is a Legendre polynomial.
+
+    The computation is done in two steps:
+
+    First, the DLT coefficients are transformed to the coefficients of a
+    discrete cosine transform (DCT) using a recursive version of the matrix
+    relation given by [1]_.
+
+    Second, the function values are computed using a DCT-III for the open
+    interval, or a DCT-I for the closed interval.
+
+    References
+    ----------
+    .. [1] Alpert, B. K., & Rokhlin, V. (1991). A fast algorithm for the
+       evaluation of Legendre expansions. SIAM Journal on Scientific and
+       Statistical Computing, 12(1), 158-179.
+
+    '''
+
+    # length n of the transform
+    if np.ndim(b) != 1:
+        raise TypeError('array must be 1d')
+    n = np.shape(b)[-1]
+
+    # type of the DCT depends on open or closed interval
+    if closed:
+        dcttype = 1
+    else:
+        dcttype = 3
+
+    # this holds the DCT coefficients
+    a = np.empty(n, dtype=float)
+
+    # these are memviews on a and b for C interop
+    cdef double[::1] a_ = a
+    cdef double[::1] b_ = b
+
+    # transform DLT coefficients to DCT coefficients using C function
+    dltdct(n, &b_[0], &a_[0])
+
+    # perform the DCT
+    xi = dct(a, type=dcttype, axis=-1, norm=None)
+
+    # done
+    return xi
+
+
+def dltmtx(n, closed=False):
+    r'''discrete Legendre transform matrix
+
+    Computes a matrix that performs the discrete Legendre transform
+    :func:`flt.dlt` when multiplied by a vector of function values.
+
+    Parameters
+    ----------
+    n : int
+        Length of the transform.
+    closed : bool, optional
+        Compute DLT over open (``closed=False``) or closed (``closed=True``)
+        interval.
+
+    Returns
+    -------
+    m : (n, n) array_like
+        Discrete Legendre transform matrix.
+
+    See Also
+    --------
+    flt.dlt : the equivalent operation
+    flt.theta : compute the angles at which the function is evaluated
+
+    Notes
+    -----
+    The discrete Legendre transform :func:`flt.dlt` performs the transformation
+    in place and does not compute the matrix :func:`flt.dltmtx`.
+
+    '''
+
+    # type of the DCT depends on open or closed interval
+    if closed:
+        dcttype = 1
+    else:
+        dcttype = 3
+
+    # compute the DCT matrix row by row (not column by column)
+    a = idct(np.eye(n), type=dcttype, axis=1, norm=None)
+
+    # this holds the DLT matrix
+    b = np.empty((n, n), dtype=float)
+
+    # these are memviews on a and b for C interop
+    cdef double[:, ::1] a_ = a
+    cdef double[:, ::1] b_ = b
+
+    # transform DCT row to DLT row using C function
+    for i in range(n):
+        dctdlt(n, &a_[i, 0], &b_[i, 0])
+
+    # return transpose since we worked on rows, not columns
+    return b.T
+
+
+def idltmtx(n, closed=False):
+    r'''inverse discrete Legendre transform matrix
+
+    Computes a matrix that performs the inverse discrete Legendre transform
+    :func:`flt.idlt` when multiplied by a vector of coefficients.
+
+    Parameters
+    ----------
+    n : int
+        Length of the transform.
+    closed : bool, optional
+        Compute inverse DLT over open (``closed=False``) or closed
+        (``closed=True``) interval.
+
+    Returns
+    -------
+    m : (n, n) array_like
+        Inverse discrete Legendre transform matrix.
+
+    See Also
+    --------
+    flt.idlt : the equivalent operation
+    flt.theta : compute the angles at which the function is evaluated
+
+    Notes
+    -----
+    The inverse discrete Legendre transform :func:`flt.idlt` performs the
+    transformation in place and does not compute the matrix
+    :func:`flt.idltmtx`.
+
+    '''
+
+    # type of the DCT depends on open or closed interval
+    if closed:
+        dcttype = 1
+    else:
+        dcttype = 3
+
+    # this is the input matrix
+    b = np.eye(n, dtype=float)
+
+    # this holds the DLT part of the matrix
+    a = np.empty((n, n), dtype=float)
+
+    # these are memviews on a and b for C interop
+    cdef double[:, ::1] a_ = a
+    cdef double[:, ::1] b_ = b
+
+    # transform DLT row to DCT row using C function
+    for i in range(n):
+        dltdct(n, &b_[i, 0], &a_[i, 0])
+
+    # multiply by DCT matrix
+    # return transpose since we worked on rows, not columns
+    return dct(a, type=dcttype, axis=1, norm=None).T
+
+
+def theta(n, closed=False):
+    r'''compute angles for DLT function values
+
+    Returns :math:`n` angles :math:`\theta_0, \ldots, \theta_{n-1}` at which
+    the function :math:`f(\cos\theta)` is evaluated in :func:`flt.dlt` or
+    :func:`flt.idlt`.
+
+    The returned angles can be distributed either over the open interval
+    :math:`(0, \theta)`, or over the closed interval :math:`[0, \pi]`, in which
+    case :math:`\theta_0 = 0, \theta_{n-1} = \pi`.
+
+    Parameters
+    ----------
+    n : int
+        Number of nodes.
+
+    Returns
+    -------
+    theta : array_like (n,)
+        Angles in radians.
+    closed : bool, optional
+        Compute angles over open (``closed=False``) or closed (``closed=True``)
+        interval.
+
+    '''
+
+    if closed:
+        t = np.linspace(0, np.pi, n, dtype=float)
+    else:
+        t = np.arange(n, dtype=float)
+        t += 0.5
+        t *= np.pi/n
+
+    return t
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 0000000..bc2be67
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,2 @@
+[build-system]
+requires = ["setuptools", "wheel", "Cython"]
diff --git a/setup.cfg b/setup.cfg
new file mode 100644
index 0000000..b99619f
--- /dev/null
+++ b/setup.cfg
@@ -0,0 +1,21 @@
+[metadata]
+name = flt
+version = 2021.3.18
+maintainer = Nicolas Tessore
+maintainer_email = n.tessore@ucl.ac.uk
+description = fast Legendre transform
+long_description = file: README.md
+long_description_content_type = text/markdown
+license = MIT
+license_file = LICENSE
+url = https://github.com/ntessore/flt
+classifiers = 
+    Programming Language :: Python :: 3
+    License :: OSI Approved :: MIT License
+    Operating System :: OS Independent
+
+[options]
+python_requires = >=3.6
+install_requires =
+    numpy
+    scipy
diff --git a/setup.py b/setup.py
new file mode 100644
index 0000000..bc0a058
--- /dev/null
+++ b/setup.py
@@ -0,0 +1,6 @@
+from setuptools import setup
+from Cython.Build import cythonize
+
+setup(
+    ext_modules=cythonize('flt.pyx'),
+)