diff --git a/bin/all_rst_to_pxd.sh b/bin/all_rst_to_pxd.sh
index 091d755e..019fb667 100755
--- a/bin/all_rst_to_pxd.sh
+++ b/bin/all_rst_to_pxd.sh
@@ -72,7 +72,7 @@ modules=(
     # "dlog"
     # "bool_mat"
     # "perm"
-    # "qfb"
+    "qfb"
     # "nf"
     # "nf_elem"
     # "fmpzi"
diff --git a/setup.py b/setup.py
index 27163c17..3780cbc8 100644
--- a/setup.py
+++ b/setup.py
@@ -120,6 +120,7 @@
     ("flint.types.acb_series", ["src/flint/types/acb_series.pyx"]),
 
     ("flint.types.dirichlet", ["src/flint/types/dirichlet.pyx"]),
+    ("flint.types.qfb", ["src/flint/types/qfb.pyx"]),
 
     ("flint.functions.showgood", ["src/flint/functions/showgood.pyx"]),
 ]
diff --git a/src/flint/__init__.py b/src/flint/__init__.py
index bf055b80..53c5bc3c 100644
--- a/src/flint/__init__.py
+++ b/src/flint/__init__.py
@@ -39,6 +39,7 @@
 from .types.acb_mat import acb_mat
 from .types.acb_series import acb_series
 
+from .types.qfb import qfb
 from .types.dirichlet import dirichlet_char, dirichlet_group
 from .functions.showgood import good, showgood
 
diff --git a/src/flint/flintlib/functions/qfb.pxd b/src/flint/flintlib/functions/qfb.pxd
new file mode 100644
index 00000000..4389e8c1
--- /dev/null
+++ b/src/flint/flintlib/functions/qfb.pxd
@@ -0,0 +1,35 @@
+from flint.flintlib.types.flint cimport fmpz_t, slong, ulong
+from flint.flintlib.types.qfb cimport qfb_t
+
+# unknown type qfb
+# unknown type qfb_hash_t
+
+
+cdef extern from "flint/qfb.h":
+    void qfb_init(qfb_t q)
+    void qfb_clear(qfb_t q)
+    # void qfb_array_clear(qfb ** forms, slong num)
+    # qfb_hash_t * qfb_hash_init(slong depth)
+    # void qfb_hash_clear(qfb_hash_t * qhash, slong depth)
+    # void qfb_hash_insert(qfb_hash_t * qhash, qfb_t q, qfb_t q2, slong iter, slong depth)
+    # slong qfb_hash_find(qfb_hash_t * qhash, qfb_t q, slong depth)
+    void qfb_set(qfb_t f, qfb_t g)
+    int qfb_equal(qfb_t f, qfb_t g)
+    void qfb_print(qfb_t q)
+    void qfb_discriminant(fmpz_t D, qfb_t f)
+    void qfb_reduce(qfb_t r, qfb_t f, fmpz_t D)
+    int qfb_is_reduced(qfb_t r)
+    # slong qfb_reduced_forms(qfb ** forms, slong d)
+    # slong qfb_reduced_forms_large(qfb ** forms, slong d)
+    void qfb_nucomp(qfb_t r, const qfb_t f, const qfb_t g, fmpz_t D, fmpz_t L)
+    void qfb_nudupl(qfb_t r, const qfb_t f, fmpz_t D, fmpz_t L)
+    void qfb_pow_ui(qfb_t r, qfb_t f, fmpz_t D, ulong exp)
+    void qfb_pow(qfb_t r, qfb_t f, fmpz_t D, fmpz_t exp)
+    void qfb_inverse(qfb_t r, qfb_t f)
+    int qfb_is_principal_form(qfb_t f, fmpz_t D)
+    void qfb_principal_form(qfb_t f, fmpz_t D)
+    int qfb_is_primitive(qfb_t f)
+    void qfb_prime_form(qfb_t r, fmpz_t D, fmpz_t p)
+    int qfb_exponent_element(fmpz_t exponent, qfb_t f, fmpz_t n, ulong B1, ulong B2_sqrt)
+    int qfb_exponent(fmpz_t exponent, fmpz_t n, ulong B1, ulong B2_sqrt, slong c)
+    int qfb_exponent_grh(fmpz_t exponent, fmpz_t n, ulong B1, ulong B2_sqrt)
diff --git a/src/flint/flintlib/types/qfb.pxd b/src/flint/flintlib/types/qfb.pxd
new file mode 100644
index 00000000..fc140b4c
--- /dev/null
+++ b/src/flint/flintlib/types/qfb.pxd
@@ -0,0 +1,10 @@
+from flint.flintlib.types.flint cimport fmpz_t
+
+cdef extern from "flint/qfb.h":
+
+    ctypedef struct qfb_struct:
+        fmpz_t a
+        fmpz_t b
+        fmpz_t c
+
+    ctypedef qfb_struct qfb_t[1]
diff --git a/src/flint/test/test_all.py b/src/flint/test/test_all.py
index cdd3f06d..714f991c 100644
--- a/src/flint/test/test_all.py
+++ b/src/flint/test/test_all.py
@@ -1646,6 +1646,41 @@ def test_nmod_series():
     pass
 
 
+def test_qfb():
+    Q = flint.qfb
+
+    assert raises(lambda: Q(1, 2, "asd"), TypeError)
+    assert raises(lambda: Q(1, "asd", 2), TypeError)
+    assert raises(lambda: Q("asd", 1, 2), TypeError)
+
+    q = Q.prime_form(-163, 53)
+    assert repr(q) == "qfb(53, 7, 1)"
+    assert q == q
+    assert q != "a"
+    assert q == Q(53, 7, 1)
+    assert not q.is_reduced()
+    assert q.reduce() == Q(1, 1, 41)
+
+    q = Q.prime_form(-199, 2)
+    assert q.is_reduced()
+    assert q**0 == Q(1, 1, 50)
+    assert q**9 == Q(1, 1, 50)
+    assert q**2 * q**5 == q**7
+    assert q.inverse() == q**-1
+    assert q.inverse() == q**8
+
+    assert raises(lambda: pow(q, 2, "asd"), NotImplementedError)
+    assert raises(lambda: q**"asd", TypeError)
+
+    q = Q.prime_form(-3212123, 7)
+    assert q**123456789123456789123456789123456789 == q.inverse()
+    assert q**-123456789123456789123456789123456789 == q
+
+    q = Q(291233996924844144901, 405366016683999883959, 141056340620716310090)
+    assert q.discriminant() == -976098765432101234567890679
+    assert q**18045470076579 == Q(1, 1, 244024691358025308641972670)
+
+
 def test_arb():
     A = flint.arb
     assert A(3) > A(2.5)
@@ -5037,6 +5072,8 @@ def test_all_tests():
     test_fq_default,
     test_fq_default_poly,
 
+    test_qfb,
+
     test_arb,
 
     test_pickling,
diff --git a/src/flint/types/meson.build b/src/flint/types/meson.build
index b5c2f953..b14547d0 100644
--- a/src/flint/types/meson.build
+++ b/src/flint/types/meson.build
@@ -30,6 +30,8 @@ exts = [
   'fq_default',
   'fq_default_poly',
 
+  'qfb',
+
   'arf',
 
   'arb',
diff --git a/src/flint/types/qfb.pxd b/src/flint/types/qfb.pxd
new file mode 100644
index 00000000..8f0df4ee
--- /dev/null
+++ b/src/flint/types/qfb.pxd
@@ -0,0 +1,7 @@
+from flint.flintlib.functions.qfb cimport *
+from flint.types.fmpz cimport fmpz
+
+cdef class qfb:
+    cdef qfb_t val
+    # the discriminant, stored for convenience
+    cdef fmpz D
diff --git a/src/flint/types/qfb.pyx b/src/flint/types/qfb.pyx
new file mode 100644
index 00000000..7f8dca7c
--- /dev/null
+++ b/src/flint/types/qfb.pyx
@@ -0,0 +1,132 @@
+from flint.flintlib.functions.fmpz cimport fmpz_abs, fmpz_root, fmpz_set
+from flint.types.fmpz cimport fmpz, any_as_fmpz
+from flint.utils.typecheck cimport typecheck
+
+cdef class qfb:
+    """
+    The qfb type represents definite binary quadratic forms
+    over Z, with composition, inverse and power operations
+    compatible with the class group of a given discriminant.
+
+    Some operations require the form to be primitive.
+    """
+    def __cinit__(self):
+        qfb_init(self.val)
+        self.D = fmpz(0)
+
+    def __dealloc__(self):
+        qfb_clear(self.val)
+
+    def __init__(self, a, b, c):
+        a_fmpz = any_as_fmpz(a)
+        b_fmpz = any_as_fmpz(b)
+        c_fmpz = any_as_fmpz(c)
+        if a_fmpz is NotImplemented:
+            raise TypeError(f"Incorrect type {type(a)} for qfb coefficient")
+        if b_fmpz is NotImplemented:
+            raise TypeError(f"Incorrect type {type(b)} for qfb coefficient")
+        if c_fmpz is NotImplemented:
+            raise TypeError(f"Incorrect type {type(c)} for qfb coefficient")
+        fmpz_set(self.val[0].a, (<fmpz>a_fmpz).val)
+        fmpz_set(self.val[0].b, (<fmpz>b_fmpz).val)
+        fmpz_set(self.val[0].c, (<fmpz>c_fmpz).val)
+        D = fmpz()
+        qfb_discriminant(D.val, self.val)
+        self.D = D
+
+    def __repr__(self):
+        a, b, c = self.coefficients()
+        return f"qfb({a}, {b}, {c})"
+
+    def __eq__(self, other):
+        if self is other:
+            return True
+
+        if typecheck(other, qfb):
+            return bool(qfb_equal(self.val, (<qfb>other).val))
+
+        return False
+
+    def __mul__(q1, q2):
+        "Returns a reduced form equivalent to the composition of q1 and q2"
+        if not q1.is_primitive():
+            raise ValueError(f"{q1} is not primitive")
+
+        cdef qfb res = qfb.__new__(qfb)
+        cdef fmpz_t L
+        fmpz_abs(L, q1.D.val)
+        fmpz_root(L, L, 4)
+        qfb_nucomp(res.val, q1.val, (<qfb>q2).val, q1.D.val, L)
+        qfb_reduce(res.val, res.val, q1.D.val)
+        res.D = q1.D
+        return res
+
+    def __pow__(q, e, mod):
+        "Returns a reduced form equivalent to the e-th iterated composition of q"
+        if mod is not None:
+            raise NotImplementedError("modular exponentiation")
+
+        if not q.is_primitive():
+            raise ValueError(f"{q} is not primitive")
+
+        e_fmpz = any_as_fmpz(e)
+        if e_fmpz is NotImplemented:
+            raise TypeError(f"exponent cannot be cast to an fmpz type: {e}")
+
+        # qfb_pow does not support negative exponents and will loop forever
+        # if a negative integer is provided.
+        e_abs = abs(e_fmpz)
+
+        cdef qfb res = qfb.__new__(qfb)
+        qfb_pow(res.val, q.val, q.D.val, (<fmpz>e_abs).val)
+        if e_fmpz < 0:
+            qfb_inverse(res.val, res.val)
+        res.D = q.D
+        return res
+
+    def coefficients(self):
+        """
+        Returns coefficients (a, b, c) of the form as a polynomial q(x,y)=ax²+bxy+cy²
+        """
+        a = fmpz()
+        fmpz_set(a.val, self.val[0].a)
+        b = fmpz()
+        fmpz_set(b.val, self.val[0].b)
+        c = fmpz()
+        fmpz_set(c.val, self.val[0].c)
+        return a, b, c
+
+    def discriminant(self):
+        return self.D
+
+    def is_reduced(self):
+        return bool(qfb_is_reduced(self.val))
+
+    def is_primitive(self):
+        return bool(qfb_is_primitive(self.val))
+
+    def inverse(self):
+        cdef qfb res = qfb.__new__(qfb)
+        qfb_inverse(res.val, self.val)
+        res.D = self.D
+        return res
+
+    def reduce(self):
+        cdef qfb res = qfb.__new__(qfb)
+        qfb_reduce(res.val, self.val, self.D.val)
+        res.D = self.D
+        return res
+
+    @classmethod
+    def prime_form(cls, D, p):
+        """
+        Returns the unique reduced form with 0 < b ≤ p. Requires that p is prime.
+        """
+
+        d_fmpz = any_as_fmpz(D)
+        p_fmpz = any_as_fmpz(p)
+
+        cdef qfb res = qfb.__new__(qfb)
+        qfb_prime_form(res.val, (<fmpz>d_fmpz).val, (<fmpz>p_fmpz).val)
+        res.D = d_fmpz
+        return res