Tests to see if ideal in quaternion algebra is primitive (cyclic)

src/doc/en/reference/references/index.rst

@@ -6342,6 +6342,8 @@ REFERENCES:
 .. [Voi2012] \J. Voight. Identifying the matrix ring: algorithms for
              quaternion algebras and quadratic forms, to appear.
 
+.. [Voi2021] \J. Voight. Quaternion algebras, Springer Nature (2021).
+
 .. [VS06]   \G.D. Villa Salvador. Topics in the Theory of Algebraic Function
             Fields. Birkh\"auser, 2006.
 

src/sage/algebras/quatalg/quaternion_algebra.py

@@ -56,7 +56,7 @@
 from sage.structure.category_object import normalize_names
 from sage.structure.parent import Parent
 from sage.matrix.matrix_space import MatrixSpace
-from sage.matrix.constructor import diagonal_matrix, matrix
+from sage.matrix.constructor import diagonal_matrix, matrix, block_matrix
 from sage.structure.sequence import Sequence
 from sage.structure.element import is_RingElement
 from sage.structure.factory import UniqueFactory
@@ -3059,6 +3059,114 @@
             ans.append(J)
         return ans
 
+    def is_integral(self):
+        r"""
+        Check if a quaternion fractional ideal is integral. An ideal in a quaternion algebra is
+        said integral if it is contained in its left order. If the left order is already defined it just
+        check the definition, otherwise it uses one of the alternative definition of Lemma 16.2.8 of
+        [Voi2021]_.
+
+        OUTPUT: a boolean.
+
+        EXAMPLES::
+
+            sage: R.<i,j,k> = QuaternionAlgebra(QQ, -1,-11)
+            sage: I = R.ideal([2 + 2*j + 140*k, 2*i + 4*j + 150*k, 8*j + 104*k, 152*k])
+            sage: I.is_integral()
+            True
+            sage: O = I.left_order()
+            sage: I.is_integral()
+            True
+            sage: I = R.ideal([1/2 + 2*j + 140*k, 2*i + 4*j + 150*k, 8*j + 104*k, 152*k])
+            sage: I.is_integral()
+            False
+
+        """
+        if self.__left_order is not None:
+            return self.free_module() <= self.left_order().free_module()
+        elif self.__right_order is not None:
+            return self.free_module() <= self.right_order().free_module()
+        else:
+            self_square = self**2
+            return self_square.free_module() <= self.free_module()
+
+    def primitive_decomposition(self):
+        r"""
+        Let `I` = ``self``. If `I` is an integral left `\mathcal{O}`-ideal return its decomposition
+        as an equivalent primitive ideal and an integer such that their product is the initial ideal.
+
+        OUTPUTS: and quivalent primitive ideal to `I`, i.e. equivalent ideal not contained in `n\mathcal{O}` for any `n>0`, and the smallest integer such that `I \subset g\mathcal{O}`.
+
+        EXAMPLES::
+
+            sage: A.<i,j,k> = QuaternionAlgebra(QQ, -1,-11)
+            sage: I = A.ideal([1/2 + 1/2*i + 1/2*j + 3/2*k, i + k, j + k, 2*k])
+            sage: I.primitive_decomposition()
+            (Fractional ideal (1/2 + 1/2*i + 1/2*j + 3/2*k, i + k, j + k, 2*k), 1)
+            sage: J = A.ideal([7/2 + 7/2*i + 49/2*j + 91/2*k, 7*i + 21*k, 35*j + 35*k, 70*k])
+            sage: Jequiv, g = J.primitive_decomposition()
+            sage: Jequiv*g == J
+            True
+            sage: Jequiv, g
+            (Fractional ideal (1/2 + 1/2*i + 7/2*j + 13/2*k, i + 3*k, 5*j + 5*k, 10*k), 7)
+
+        TESTS::
+
+        Checks on random crafted ideals that they decompose as expected.
+
+            sage: for d in ( m for m in range(400, 750) if is_squarefree(m) ):
+            ....:     A = QuaternionAlgebra(d)
+            ....:     O = A.maximal_order()
+            ....:     for _ in range(10):
+            ....:         a = O.random_element()
+            ....:         if not a.is_constant(): # avoids a = 0
+            ....:             I = a*O + a.reduced_norm()*O
+            ....:             if I.is_integral():
+            ....:                 J,g = I.primitive_decomposition()
+            ....:                 assert J*g == I
+            ....:                 assert J.is_primitive()
+
+        """
+
+        if not self.is_integral():
+            raise ValueError("primitive ideals are defined only for integral ideals")
+
+        I_basis = self.basis_matrix()
+        O_basis = self.left_order().basis_matrix()
+
+        # Write I in the basis of its left order via rref
+        M = O_basis.solve_left(I_basis)
+        g = Integer(gcd(M.list()))
+
+        # If g is 1 then the ideal is primitive
+        if g.is_one():
+            return self, g
+
+        J = self.scale(1/g)
+
+        return J, g
+
+    def is_primitive(self):
+        r"""
+        Check if the quaternion fractional ideal is primitive. An integral left
+        $O$-ideal for some order $O$ is said primitive if for all integers $n > 1$
+        $I$ is not contained in $nO$.
+
+        OUTPUT: a boolean.
+
+        EXAMPLES::
+
+            sage: A.<i,j,k> = QuaternionAlgebra(QQ, -1,-11)
+            sage: I = A.ideal([1/2 + 1/2*i + 1/2*j + 3/2*k, i + k, j + k, 2*k])
+            sage: I.is_primitive()
+            True
+            sage: (2*I).is_primitive()
+            False
+
+        """
+        _,g = self.primitive_decomposition()
+        return g.is_one()
+
 #######################################################################
 # Some utility functions that are needed here and are too
 # specialized to go elsewhere.