Add bindings for type qfb_t

bin/all_rst_to_pxd.sh

@@ -72,7 +72,7 @@ modules=(
     # "dlog"
     # "bool_mat"
     # "perm"
-    # "qfb"
+    "qfb"
     # "nf"
     # "nf_elem"
     # "fmpzi"

src/flint/__init__.py

@@ -26,6 +26,8 @@
 
 from .types.fmpq_mpoly import fmpq_mpoly_ctx, fmpq_mpoly, fmpq_mpoly_vec
 
+from .types.qfb import *
+
 from .types.fq_default import *
 from .types.fq_default_poly import *
 

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)

src/flint/flintlib/types/qfb.pxd

@@ -0,0 +1,11 @@
+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]
+

src/flint/test/test_all.py

@@ -1617,6 +1617,34 @@ def test_nmod_series():
     # XXX: currently no code in nmod_series.pyx
     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(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
+
+    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
@@ -4897,6 +4925,8 @@ def test_all_tests():
     test_fq_default,
     test_fq_default_poly,
 
+    test_qfb,
+
     test_arb,
 
     test_pickling,

src/flint/types/meson.build

@@ -30,6 +30,8 @@ exts = [
   'fq_default',
   'fq_default_poly',
 
+  'qfb',
+
   'arf',
 
   'arb',

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

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