text: """
  Quaternion Algebras
  
  AUTHORS:
  
  - Jon Bobber (2009): rewrite
  
  - William Stein (2009): rewrite
  
  - Julian Rueth (2014-03-02): use UniqueFactory for caching
  
  - Peter Bruin (2021): do not require the base ring to be a field
  
  This code is partly based on Sage code by David Kohel from 2005.
  
  TESTS:
  
  Pickling test::
  
      sage: Q.<i,j,k> = QuaternionAlgebra(QQ,-5,-2)
      sage: Q == loads(dumps(Q))
      True
  """
  
  # ****************************************************************************
  #       Copyright (C) 2009 William Stein <wstein@gmail.com>
  #       Copyright (C) 2009 Jonathan Bober <jwbober@gmail.com>
  #       Copyright (C) 2014 Julian Rueth <julian.rueth@fsfe.org>
  #       Copyright (C) 2021 Peter Bruin <P.J.Bruin@math.leidenuniv.nl>
  #
  # This program is free software: you can redistribute it and/or modify
  # it under the terms of the GNU General Public License as published by
  # the Free Software Foundation, either version 2 of the License, or
  # (at your option) any later version.
  #                  https://www.gnu.org/licenses/
  # ****************************************************************************
  
  from sage.arith.all import (hilbert_conductor_inverse, hilbert_conductor,
                              factor, gcd, kronecker_symbol, valuation)
  from sage.rings.all import RR, Integer
  from sage.rings.integer_ring import ZZ
  from sage.rings.rational import Rational
  from sage.rings.finite_rings.finite_field_constructor import GF
  
  from sage.rings.ring import Algebra
  from sage.rings.ideal import Ideal_fractional
  from sage.rings.rational_field import is_RationalField, QQ
  from sage.rings.infinity import infinity
  from sage.rings.number_field.number_field import is_NumberField
  from sage.rings.power_series_ring import PowerSeriesRing
  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.structure.sequence import Sequence
  from sage.structure.element import is_RingElement
  from sage.structure.factory import UniqueFactory
  from sage.modules.free_module import FreeModule
  from sage.modules.free_module_element import vector
  
  from operator import itemgetter
  
  from .quaternion_algebra_element import (
      QuaternionAlgebraElement_abstract,
      QuaternionAlgebraElement_generic,
      QuaternionAlgebraElement_rational_field,
      QuaternionAlgebraElement_number_field)
  from . import quaternion_algebra_cython
  
  from sage.modular.modsym.p1list import P1List
  
  from sage.misc.cachefunc import cached_method
  
  from sage.categories.algebras import Algebras
  
  ########################################################
  # Constructor
  ########################################################
  
  
  class QuaternionAlgebraFactory(UniqueFactory):
      r"""
      Construct a quaternion algebra.
  
      INPUT:
  
      There are three input formats:
  
      - ``QuaternionAlgebra(a, b)``, where `a` and `b` can be coerced to
        units in a common field `K` of characteristic different from 2.
  
      - ``QuaternionAlgebra(K, a, b)``, where `K` is a ring in which 2
        is a unit and `a` and `b` are units of `K`.
  
      - ``QuaternionAlgebra(D)``, where `D \ge 1` is a squarefree
        integer.  This constructs a quaternion algebra of discriminant
        `D` over `K = \QQ`.  Suitable nonzero rational numbers `a`, `b`
        as above are deduced from `D`.
  
      OUTPUT:
  
      The quaternion algebra `(a, b)_K` over `K` generated by `i`, `j`
      subject to `i^2 = a`, `j^2 = b`, and `ji = -ij`.
  
      EXAMPLES:
  
      ``QuaternionAlgebra(a, b)`` -- return the quaternion algebra
      `(a, b)_K`, where the base ring `K` is a suitably chosen field
      containing `a` and `b`::
  
          sage: QuaternionAlgebra(-2,-3)
          Quaternion Algebra (-2, -3) with base ring Rational Field
          sage: QuaternionAlgebra(GF(5)(2), GF(5)(3))
          Quaternion Algebra (2, 3) with base ring Finite Field of size 5
          sage: QuaternionAlgebra(2, GF(5)(3))
          Quaternion Algebra (2, 3) with base ring Finite Field of size 5
          sage: QuaternionAlgebra(QQ[sqrt(2)](-1), -5)
          Quaternion Algebra (-1, -5) with base ring Number Field in sqrt2 with defining polynomial x^2 - 2 with sqrt2 = 1.414213562373095?
          sage: QuaternionAlgebra(sqrt(-1), sqrt(-3))
          Quaternion Algebra (I, sqrt(-3)) with base ring Symbolic Ring
          sage: QuaternionAlgebra(1r,1)
          Quaternion Algebra (1, 1) with base ring Rational Field
          sage: A.<t> = ZZ[]
          sage: QuaternionAlgebra(-1, t)
          Quaternion Algebra (-1, t) with base ring Fraction Field of Univariate Polynomial Ring in t over Integer Ring
  
      Python ints and floats may be passed to the
      ``QuaternionAlgebra(a, b)`` constructor, as may all pairs of
      nonzero elements of a domain not of characteristic 2.
  
      The following tests address the issues raised in :trac:`10601`::
  
          sage: QuaternionAlgebra(1r,1)
          Quaternion Algebra (1, 1) with base ring Rational Field
          sage: QuaternionAlgebra(1,1.0r)
          Quaternion Algebra (1.00000000000000, 1.00000000000000) with base ring Real Field with 53 bits of precision
          sage: QuaternionAlgebra(0,0)
          Traceback (most recent call last):
          ...
          ValueError: defining elements of quaternion algebra (0, 0) are not invertible in Rational Field
          sage: QuaternionAlgebra(GF(2)(1),1)
          Traceback (most recent call last):
          ...
          ValueError: 2 is not invertible in Finite Field of size 2
          sage: a = PermutationGroupElement([1,2,3])
          sage: QuaternionAlgebra(a, a)
          Traceback (most recent call last):
          ...
          ValueError: a and b must be elements of a ring with characteristic not 2
  
      ``QuaternionAlgebra(K, a, b)`` -- return the quaternion algebra
      defined by `(a, b)` over the ring `K`::
  
          sage: QuaternionAlgebra(QQ, -7, -21)
          Quaternion Algebra (-7, -21) with base ring Rational Field
          sage: QuaternionAlgebra(QQ[sqrt(2)], -2,-3)
          Quaternion Algebra (-2, -3) with base ring Number Field in sqrt2 with defining polynomial x^2 - 2 with sqrt2 = 1.414213562373095?
  
      ``QuaternionAlgebra(D)`` -- `D` is a squarefree integer; return a
      rational quaternion algebra of discriminant `D`::
  
          sage: QuaternionAlgebra(1)
          Quaternion Algebra (-1, 1) with base ring Rational Field
          sage: QuaternionAlgebra(2)
          Quaternion Algebra (-1, -1) with base ring Rational Field
          sage: QuaternionAlgebra(7)
          Quaternion Algebra (-1, -7) with base ring Rational Field
          sage: QuaternionAlgebra(2*3*5*7)
          Quaternion Algebra (-22, 210) with base ring Rational Field
  
      If the coefficients `a` and `b` in the definition of the quaternion
      algebra are not integral, then a slower generic type is used for
      arithmetic::
  
          sage: type(QuaternionAlgebra(-1,-3).0)
          <... 'sage.algebras.quatalg.quaternion_algebra_element.QuaternionAlgebraElement_rational_field'>
          sage: type(QuaternionAlgebra(-1,-3/2).0)
          <... 'sage.algebras.quatalg.quaternion_algebra_element.QuaternionAlgebraElement_generic'>
  
      Make sure caching is sane::
  
          sage: A = QuaternionAlgebra(2,3); A
          Quaternion Algebra (2, 3) with base ring Rational Field
          sage: B = QuaternionAlgebra(GF(5)(2),GF(5)(3)); B
          Quaternion Algebra (2, 3) with base ring Finite Field of size 5
          sage: A is QuaternionAlgebra(2,3)
          True
          sage: B is QuaternionAlgebra(GF(5)(2),GF(5)(3))
          True
          sage: Q = QuaternionAlgebra(2); Q
          Quaternion Algebra (-1, -1) with base ring Rational Field
          sage: Q is QuaternionAlgebra(QQ,-1,-1)
          True
          sage: Q is QuaternionAlgebra(-1,-1)
          True
          sage: Q.<ii,jj,kk> = QuaternionAlgebra(15); Q.variable_names()
          ('ii', 'jj', 'kk')
          sage: QuaternionAlgebra(15).variable_names()
          ('i', 'j', 'k')
  
      TESTS:
  
      Verify that bug found when working on :trac:`12006` involving coercing
      invariants into the base field is fixed::
  
          sage: Q = QuaternionAlgebra(-1,-1); Q
          Quaternion Algebra (-1, -1) with base ring Rational Field
          sage: parent(Q._a)
          Rational Field
          sage: parent(Q._b)
          Rational Field
      """
      def create_key(self, arg0, arg1=None, arg2=None, names='i,j,k'):
          """
          Create a key that uniquely determines a quaternion algebra.
  
          TESTS::
  
              sage: QuaternionAlgebra.create_key(-1,-1)
              (Rational Field, -1, -1, ('i', 'j', 'k'))
          """
          # QuaternionAlgebra(D)
          if arg1 is None and arg2 is None:
              K = QQ
              D = Integer(arg0)
              a, b = hilbert_conductor_inverse(D)
              a = Rational(a)
              b = Rational(b)
  
          elif arg2 is None:
              # If arg0 or arg1 are Python data types, coerce them
              # to the relevant Sage types. This is a bit inelegant.
              L = []
              for a in [arg0, arg1]:
                  if is_RingElement(a):
                      L.append(a)
                  elif isinstance(a, int):
                      L.append(Integer(a))
                  elif isinstance(a, float):
                      L.append(RR(a))
                  else:
                      raise ValueError("a and b must be elements of a ring with characteristic not 2")
  
              # QuaternionAlgebra(a, b)
              v = Sequence(L)
              K = v.universe().fraction_field()
              a = K(v[0])
              b = K(v[1])
  
          # QuaternionAlgebra(K, a, b)
          else:
              K = arg0
              a = K(arg1)
              b = K(arg2)
  
          if not K(2).is_unit():
              raise ValueError("2 is not invertible in %s" % K)
          if not (a.is_unit() and b.is_unit()):
              raise ValueError("defining elements of quaternion algebra (%s, %s) are not invertible in %s"
                               % (a, b, K))
  
          names = normalize_names(3, names)
          return (K, a, b, names)
  
      def create_object(self, version, key, **extra_args):
          """
          Create the object from the key (extra arguments are ignored). This is
          only called if the object was not found in the cache.
  
          TESTS::
  
              sage: QuaternionAlgebra.create_object("6.0", (QQ, -1, -1, ('i', 'j', 'k')))
              Quaternion Algebra (-1, -1) with base ring Rational Field
  
          """
          K, a, b, names = key
          return QuaternionAlgebra_ab(K, a, b, names=names)
  
  
  QuaternionAlgebra = QuaternionAlgebraFactory("QuaternionAlgebra")
  
  ########################################################
  # Classes
  ########################################################
  
  
  def is_QuaternionAlgebra(A):
      """
      Return ``True`` if ``A`` is of the QuaternionAlgebra data type.
  
      EXAMPLES::
  
          sage: sage.algebras.quatalg.quaternion_algebra.is_QuaternionAlgebra(QuaternionAlgebra(QQ,-1,-1))
          True
          sage: sage.algebras.quatalg.quaternion_algebra.is_QuaternionAlgebra(ZZ)
          False
      """
      return isinstance(A, QuaternionAlgebra_abstract)
  
  
  class QuaternionAlgebra_abstract(Algebra):
      def _repr_(self):
          """
          EXAMPLES::
  
              sage: sage.algebras.quatalg.quaternion_algebra.QuaternionAlgebra_abstract(QQ)._repr_()
              'Quaternion Algebra with base ring Rational Field'
          """
          return "Quaternion Algebra with base ring %s" % self.base_ring()
  
      def ngens(self):
          """
          Return the number of generators of the quaternion algebra as a K-vector
          space, not including 1.
  
          This value is always 3: the algebra is spanned
          by the standard basis `1`, `i`, `j`, `k`.
  
          EXAMPLES::
  
              sage: Q.<i,j,k> = QuaternionAlgebra(QQ,-5,-2)
              sage: Q.ngens()
              3
              sage: Q.gens()
              [i, j, k]
          """
          return 3
  
      @cached_method
      def basis(self):
          """
          Return the fixed basis of ``self``, which is `1`, `i`, `j`, `k`, where
          `i`, `j`, `k` are the generators of ``self``.
  
          EXAMPLES::
  
              sage: Q.<i,j,k> = QuaternionAlgebra(QQ,-5,-2)
              sage: Q.basis()
              (1, i, j, k)
  
              sage: Q.<xyz,abc,theta> = QuaternionAlgebra(GF(9,'a'),-5,-2)
              sage: Q.basis()
              (1, xyz, abc, theta)
  
          The basis is cached::
  
              sage: Q.basis() is Q.basis()
              True
          """
          i, j, k = self.gens()
          return (self.one(), i, j, k)
  
      @cached_method
      def inner_product_matrix(self):
          """
          Return the inner product matrix associated to ``self``.
  
          This is the
          Gram matrix of the reduced norm as a quadratic form on ``self``.
          The standard basis `1`, `i`, `j`, `k` is orthogonal, so this matrix
          is just the diagonal matrix with diagonal entries `2`, `2a`, `2b`,
          `2ab`.
  
          EXAMPLES::
  
              sage: Q.<i,j,k> = QuaternionAlgebra(-5,-19)
              sage: Q.inner_product_matrix()
              [  2   0   0   0]
              [  0  10   0   0]
              [  0   0  38   0]
              [  0   0   0 190]
          """
          a, b = self._a, self._b
          M = diagonal_matrix(self.base_ring(), [2, -2 * a, -2 * b, 2 * a * b])
          M.set_immutable()
          return M
  
      def is_commutative(self) -> bool:
          """
          Return ``False`` always, since all quaternion algebras are
          noncommutative.
  
          EXAMPLES::
  
              sage: Q.<i,j,k> = QuaternionAlgebra(QQ, -3,-7)
              sage: Q.is_commutative()
              False
          """
          return False
  
      def is_division_algebra(self) -> bool:
          """
          Return ``True`` if the quaternion algebra is a division algebra (i.e.
          every nonzero element in ``self`` is invertible), and ``False`` if the
          quaternion algebra is isomorphic to the 2x2 matrix algebra.
  
          EXAMPLES::
  
              sage: QuaternionAlgebra(QQ,-5,-2).is_division_algebra()
              True
              sage: QuaternionAlgebra(1).is_division_algebra()
              False
              sage: QuaternionAlgebra(2,9).is_division_algebra()
              False
              sage: QuaternionAlgebra(RR(2.),1).is_division_algebra()
              Traceback (most recent call last):
              ...
              NotImplementedError: base field must be rational numbers
          """
          if not is_RationalField(self.base_ring()):
              raise NotImplementedError("base field must be rational numbers")
          return self.discriminant() != 1
  
      def is_matrix_ring(self) -> bool:
          """
          Return ``True`` if the quaternion algebra is isomorphic to the 2x2
          matrix ring, and ``False`` if ``self`` is a division algebra (i.e.
          every nonzero element in ``self`` is invertible).
  
          EXAMPLES::
  
              sage: QuaternionAlgebra(QQ,-5,-2).is_matrix_ring()
              False
              sage: QuaternionAlgebra(1).is_matrix_ring()
              True
              sage: QuaternionAlgebra(2,9).is_matrix_ring()
              True
              sage: QuaternionAlgebra(RR(2.),1).is_matrix_ring()
              Traceback (most recent call last):
              ...
              NotImplementedError: base field must be rational numbers
  
          """
          if not is_RationalField(self.base_ring()):
              raise NotImplementedError("base field must be rational numbers")
          return self.discriminant() == 1
  
      def is_exact(self) -> bool:
          """
          Return ``True`` if elements of this quaternion algebra are represented
          exactly, i.e. there is no precision loss when doing arithmetic. A
          quaternion algebra is exact if and only if its base field is
          exact.
  
          EXAMPLES::
  
              sage: Q.<i,j,k> = QuaternionAlgebra(QQ, -3, -7)
              sage: Q.is_exact()
              True
              sage: Q.<i,j,k> = QuaternionAlgebra(Qp(7), -3, -7)
              sage: Q.is_exact()
              False
          """
          return self.base_ring().is_exact()
  
      def is_field(self, proof=True) -> bool:
          """
          Return ``False`` always, since all quaternion algebras are
          noncommutative and all fields are commutative.
  
          EXAMPLES::
  
              sage: Q.<i,j,k> = QuaternionAlgebra(QQ, -3, -7)
              sage: Q.is_field()
              False
          """
          return False
  
      def is_finite(self) -> bool:
          """
          Return ``True`` if the quaternion algebra is finite as a set.
  
          Algorithm: A quaternion algebra is finite if and only if the
          base field is finite.
  
          EXAMPLES::
  
              sage: Q.<i,j,k> = QuaternionAlgebra(QQ, -3, -7)
              sage: Q.is_finite()
              False
              sage: Q.<i,j,k> = QuaternionAlgebra(GF(5), -3, -7)
              sage: Q.is_finite()
              True
          """
          return self.base_ring().is_finite()
  
      def is_integral_domain(self, proof=True) -> bool:
          """
          Return ``False`` always, since all quaternion algebras are
          noncommutative and integral domains are commutative (in Sage).
  
          EXAMPLES::
  
              sage: Q.<i,j,k> = QuaternionAlgebra(QQ, -3, -7)
              sage: Q.is_integral_domain()
              False
          """
          return False
  
      def is_noetherian(self) -> bool:
          """
          Return ``True`` always, since any quaternion algebra is a noetherian
          ring (because it is a finitely generated module over a field).
  
          EXAMPLES::
  
              sage: Q.<i,j,k> = QuaternionAlgebra(QQ, -3, -7)
              sage: Q.is_noetherian()
              True
          """
          return True
  
      def order(self):
          """
          Return the number of elements of the quaternion algebra, or
          ``+Infinity`` if the algebra is not finite.
  
          EXAMPLES::
  
              sage: Q.<i,j,k> = QuaternionAlgebra(QQ, -3, -7)
              sage: Q.order()
              +Infinity
              sage: Q.<i,j,k> = QuaternionAlgebra(GF(5), -3, -7)
              sage: Q.order()
              625
          """
          return (self.base_ring().order())**4
  
      def random_element(self, *args, **kwds):
          """
          Return a random element of this quaternion algebra.
  
          The ``args`` and ``kwds`` are passed to the ``random_element`` method
          of the base ring.
  
          EXAMPLES::
  
              sage: g = QuaternionAlgebra(QQ[sqrt(2)], -3, 7).random_element()
              sage: g.parent() is QuaternionAlgebra(QQ[sqrt(2)], -3, 7)
              True
              sage: g = QuaternionAlgebra(-3, 19).random_element()
              sage: g.parent() is QuaternionAlgebra(-3, 19)
              True
              sage: g = QuaternionAlgebra(GF(17)(2), 3).random_element()
              sage: g.parent() is QuaternionAlgebra(GF(17)(2), 3)
              True
  
          Specify the numerator and denominator bounds::
  
              sage: g = QuaternionAlgebra(-3,19).random_element(10^6, 10^6)
              sage: for h in g:
              ....:     assert h.numerator() in range(-10^6, 10^6 + 1)
              ....:     assert h.denominator() in range(10^6 + 1)
  
              sage: g = QuaternionAlgebra(-3,19).random_element(5, 4)
              sage: for h in g:
              ....:     assert h.numerator() in range(-5, 5 + 1)
              ....:     assert h.denominator() in range(4 + 1)
          """
          K = self.base_ring()
          return self([K.random_element(*args, **kwds) for _ in range(4)])
  
      @cached_method
      def free_module(self):
          """
          Return the free module associated to ``self`` with inner
          product given by the reduced norm.
  
          EXAMPLES::
  
              sage: A.<t> = LaurentPolynomialRing(GF(3))
              sage: B = QuaternionAlgebra(A, -1, t)
              sage: B.free_module()
              Ambient free quadratic module of rank 4 over the principal ideal domain Univariate Laurent Polynomial Ring in t over Finite Field of size 3
              Inner product matrix:
              [2 0 0 0]
              [0 2 0 0]
              [0 0 t 0]
              [0 0 0 t]
          """
          return FreeModule(self.base_ring(), 4, inner_product_matrix=self.inner_product_matrix())
  
      def vector_space(self):
          """
          Alias for :meth:`free_module`.
  
          EXAMPLES::
  
              sage: QuaternionAlgebra(-3,19).vector_space()
              Ambient quadratic space of dimension 4 over Rational Field
              Inner product matrix:
              [   2    0    0    0]
              [   0    6    0    0]
              [   0    0  -38    0]
              [   0    0    0 -114]
          """
          return self.free_module()
  
  
  class QuaternionAlgebra_ab(QuaternionAlgebra_abstract):
      """
      A quaternion algebra of the form `(a, b)_K`.
  
      See ``QuaternionAlgebra`` for many more examples.
  
      INPUT:
  
      - ``base_ring`` -- a commutative ring `K` in which 2 is invertible
      - ``a, b`` -- units of `K`
      - ``names`` -- string (optional, default 'i,j,k') names of the generators
  
      OUTPUT:
  
      The quaternion algebra `(a, b)` over `K` generated by `i` and `j`
      subject to `i^2 = a`, `j^2 = b`, and `ji = -ij`.
  
      EXAMPLES::
  
          sage: QuaternionAlgebra(QQ, -7, -21)  # indirect doctest
          Quaternion Algebra (-7, -21) with base ring Rational Field
      """
      def __init__(self, base_ring, a, b, names='i,j,k'):
          """
          Create the quaternion algebra with `i^2 = a`, `j^2 = b`, and
          `ij = -ji = k`.
  
          TESTS:
  
          Test making quaternion elements (using the element constructor)::
  
              sage: Q.<i,j,k> = QuaternionAlgebra(QQ,-1,-2)
              sage: a = Q(2/3); a
              2/3
              sage: type(a)
              <... 'sage.algebras.quatalg.quaternion_algebra_element.QuaternionAlgebraElement_rational_field'>
              sage: Q(a)
              2/3
              sage: Q([1,2,3,4])
              1 + 2*i + 3*j + 4*k
              sage: Q((1,2,3,4))
              1 + 2*i + 3*j + 4*k
              sage: Q(-3/5)
              -3/5
  
              sage: TestSuite(Q).run()
  
          The element 2 must be a unit in the base ring::
  
              sage: Q.<ii,jj,kk> = QuaternionAlgebra(ZZ,-5,-19)
              Traceback (most recent call last):
              ...
              ValueError: 2 is not invertible in Integer Ring
          """
          cat = Algebras(base_ring).Division().FiniteDimensional()
          Algebra.__init__(self, base_ring, names, category=cat)
          self._a = a
          self._b = b
          if is_RationalField(base_ring) and a.denominator() == 1 == b.denominator():
              self.Element = QuaternionAlgebraElement_rational_field
          elif (is_NumberField(base_ring) and base_ring.degree() > 2 and base_ring.is_absolute() and
                a.denominator() == 1 == b.denominator() and base_ring.defining_polynomial().is_monic()):
              # This QuaternionAlgebraElement_number_field class is not
              # designed to work with elements of a quadratic field.  To
              # do that, the main thing would be to implement
              # __getitem__, etc.  This would maybe give a factor of 2
              # (or more?) speedup.  Much care must be taken because the
              # underlying representation of quadratic fields is a bit
              # tricky.
              self.Element = QuaternionAlgebraElement_number_field
          else:
              self.Element = QuaternionAlgebraElement_generic
          self._populate_coercion_lists_(coerce_list=[base_ring])
          self._gens = [self([0, 1, 0, 0]), self([0, 0, 1, 0]), self([0, 0, 0, 1])]
  
      @cached_method
      def maximal_order(self, take_shortcuts=True):
          r"""
          Return a maximal order in this quaternion algebra.
  
          The algorithm used is from [Voi2012]_.
  
          INPUT:
  
          - ``take_shortcuts`` -- (default: ``True``) if the discriminant is
            prime and the invariants of the algebra are of a nice form, use
            Proposition 5.2 of [Piz1980]_.
  
          OUTPUT:
  
          A maximal order in this quaternion algebra.
  
          EXAMPLES::
  
              sage: QuaternionAlgebra(-1,-7).maximal_order()
              Order of Quaternion Algebra (-1, -7) with base ring Rational Field with basis (1/2 + 1/2*j, 1/2*i + 1/2*k, j, k)
  
              sage: QuaternionAlgebra(-1,-1).maximal_order().basis()
              (1/2 + 1/2*i + 1/2*j + 1/2*k, i, j, k)
  
              sage: QuaternionAlgebra(-1,-11).maximal_order().basis()
              (1/2 + 1/2*j, 1/2*i + 1/2*k, j, k)
  
              sage: QuaternionAlgebra(-1,-3).maximal_order().basis()
              (1/2 + 1/2*j, 1/2*i + 1/2*k, j, k)
  
              sage: QuaternionAlgebra(-3,-1).maximal_order().basis()
              (1/2 + 1/2*i, 1/2*j - 1/2*k, i, -k)
  
              sage: QuaternionAlgebra(-2,-5).maximal_order().basis()
              (1/2 + 1/2*j + 1/2*k, 1/4*i + 1/2*j + 1/4*k, j, k)
  
              sage: QuaternionAlgebra(-5,-2).maximal_order().basis()
              (1/2 + 1/2*i - 1/2*k, 1/2*i + 1/4*j - 1/4*k, i, -k)
  
              sage: QuaternionAlgebra(-17,-3).maximal_order().basis()
              (1/2 + 1/2*j, 1/2*i + 1/2*k, -1/3*j - 1/3*k, k)
  
              sage: QuaternionAlgebra(-3,-17).maximal_order().basis()
              (1/2 + 1/2*i, 1/2*j - 1/2*k, -1/3*i + 1/3*k, -k)
  
              sage: QuaternionAlgebra(-17*9,-3).maximal_order().basis()
              (1, 1/3*i, 1/6*i + 1/2*j, 1/2 + 1/3*j + 1/18*k)
  
              sage: QuaternionAlgebra(-2, -389).maximal_order().basis()
              (1/2 + 1/2*j + 1/2*k, 1/4*i + 1/2*j + 1/4*k, j, k)
  
          If you want bases containing 1, switch off ``take_shortcuts``::
  
              sage: QuaternionAlgebra(-3,-89).maximal_order(take_shortcuts=False)
              Order of Quaternion Algebra (-3, -89) with base ring Rational Field with basis (1, 1/2 + 1/2*i, j, 1/2 + 1/6*i + 1/2*j + 1/6*k)
  
              sage: QuaternionAlgebra(1,1).maximal_order(take_shortcuts=False)    # Matrix ring
              Order of Quaternion Algebra (1, 1) with base ring Rational Field with basis (1, 1/2 + 1/2*i, j, 1/2*j + 1/2*k)
  
              sage: QuaternionAlgebra(-22,210).maximal_order(take_shortcuts=False)
              Order of Quaternion Algebra (-22, 210) with base ring Rational Field with basis (1, i, 1/2*i + 1/2*j, 1/2 + 17/22*i + 1/44*k)
  
              sage: for d in ( m for m in range(1, 750) if is_squarefree(m) ):        # long time (3s)
              ....:     A = QuaternionAlgebra(d)
              ....:     R = A.maximal_order(take_shortcuts=False)
              ....:     assert A.discriminant() == R.discriminant()
  
          We do not support number fields other than the rationals yet::
  
              sage: K = QuadraticField(5)
              sage: QuaternionAlgebra(K,-1,-1).maximal_order()
              Traceback (most recent call last):
              ...
              NotImplementedError: maximal order only implemented for rational quaternion algebras
          """
          if self.base_ring() != QQ:
              raise NotImplementedError("maximal order only implemented for rational quaternion algebras")
  
          d_A = self.discriminant()
  
          # The following only works over QQ if the discriminant is prime
          # and if the invariants are of the special form
          # (every quaternion algebra of prime discriminant has a representation
          #  of such a form though)
          a, b = self.invariants()
          if take_shortcuts and d_A.is_prime() and a in ZZ and b in ZZ:
              a = ZZ(a)
              b = ZZ(b)
              i, j, k = self.gens()
  
              # if necessary, try to swap invariants to match Pizer's paper
              if (a != -1 and b == -1) or (b == -2) \
                 or (a != -1 and a != -2 and (-a) % 8 != 1):
                  a, b = b, a
                  i, j = j, i
                  k = i * j
  
              basis = []
              if (a, b) == (-1, -1):
                  basis = [(1+i+j+k)/2, i, j, k]
              elif a == -1 and (-b).is_prime() and ((-b) % 4 == 3):
                  basis = [(1+j)/2, (i+k)/2, j, k]
              elif a == -2 and (-b).is_prime() and ((-b) % 8 == 5):
                  basis = [(1+j+k)/2, (i+2*j+k)/4, j, k]
              elif (-a).is_prime() and (-b).is_prime():
                  q = -b
                  p = -a
  
                  if q % 4 == 3 and kronecker_symbol(p, q) == -1:
                      a = 0
                      while (a * a * p + 1) % q:
                          a += 1
                      basis = [(1+j)/2, (i+k)/2, -(j+a*k)/q, k]
  
              if basis:
                  return self.quaternion_order(basis)
  
          # The following code should always work (over QQ)
          # Start with <1,i,j,k>
          R = self.quaternion_order([1] + self.gens())
          d_R = R.discriminant()
  
          e_new_gens = []
  
          # For each prime at which R is not yet maximal, make it bigger
          for p, _ in d_R.factor():
              e = R.basis()
              while self.quaternion_order(e).discriminant().valuation(p) > d_A.valuation(p):
                  # Compute a normalized basis at p
                  f = normalize_basis_at_p(list(e), p)
  
                  # Ensure the basis lies in R by clearing denominators
                  # (this may make the order smaller at q != p)
                  # Also saturate the basis (divide out p as far as possible)
                  V = self.base_ring()**4
                  A = matrix(self.base_ring(), 4, 4, [list(g) for g in e])
  
                  e_n = []
                  x_rows = A.solve_left(matrix([V(vec.coefficient_tuple())
                                                for (vec, val) in f]),
                                        check=False).rows()
                  denoms = [x.denominator() for x in x_rows]
                  for i in range(4):
                      vec = f[i][0]
                      val = f[i][1]
  
                      v = (val/2).floor()
                      e_n.append(denoms[i] / p**(v) * vec)
  
                  # for e_n to become p-saturated we still need to sort by
                  # ascending valuation of the quadratic form
                  lst = sorted(zip(e_n, [f[m][1].mod(2) for m in range(4)]),
                               key=itemgetter(1))
                  e_n = list(next(zip(*lst)))
  
                  # Final step: Enlarge the basis at p
                  if p != 2:
                      # ensure that v_p(e_n[1]**2) = 0 by swapping basis elements
                      if ZZ(e_n[1]**2).valuation(p) != 0:
                          if ZZ(e_n[2]**2).valuation(p) == 0:
                              e_n[1], e_n[2] = e_n[2], e_n[1]
                          else:
                              e_n[1], e_n[3] = e_n[3], e_n[1]
  
                      a = ZZ(e_n[1]**2)
                      b = ZZ(e_n[2]**2)
  
                      if b.valuation(p) > 0:      # if v_p(b) = 0, then already p-maximal
                          F = ZZ.quo(p)
                          if F(a).is_square():
                              x = F(a).sqrt().lift()
                              if (x**2 - a).mod(p**2) == 0:  # make sure v_p(x**2 - a) = 1
                                  x = x + p
                              g = 1/p*(x - e_n[1])*e_n[2]
                              e_n[2] = g
                              e_n[3] = e_n[1]*g
  
                  else:   # p == 2
                      t = e_n[1].reduced_trace()
                      a = -e_n[1].reduced_norm()
                      b = ZZ(e_n[2]**2)
  
                      if t.valuation(p) == 0:
                          if b.valuation(p) > 0:
                              x = a
                              if (x**2 - t*x + a).mod(p**2) == 0:  # make sure v_p(...) = 1
                                  x = x + p
                              g = 1/p*(x - e_n[1])*e_n[2]
                              e_n[2] = g
                              e_n[3] = e_n[1]*g
  
                      else:   # t.valuation(p) > 0
                          (y, z, w) = maxord_solve_aux_eq(a, b, p)
                          g = 1/p*(1 + y*e_n[1] + z*e_n[2] + w*e_n[1]*e_n[2])
                          h = (z*b)*e_n[1] - (y*a)*e_n[2]
                          e_n[1:4] = [g, h, g * h]
                          if (1 - a*y**2 - b*z**2 + a*b*w**2).valuation(2) > 2:
                              e_n = basis_for_quaternion_lattice(list(e) + e_n[1:], reverse=True)
  
                  # e_n now contains elements that locally at p give a bigger order,
                  # but the basis may be messed up at other primes (it might not even
                  # be an order). We will join them all together at the end
                  e = e_n
  
              e_new_gens.extend(e[1:])
  
          e_new = basis_for_quaternion_lattice(list(R.basis()) + e_new_gens, reverse=True)
          return self.quaternion_order(e_new)
  
      def invariants(self):
          """
          Return the structural invariants `a`, `b` of this quaternion
          algebra: ``self`` is generated by `i`, `j` subject to
          `i^2 = a`, `j^2 = b` and `ji = -ij`.
  
          EXAMPLES::
  
              sage: Q.<i,j,k> = QuaternionAlgebra(15)
              sage: Q.invariants()
              (-3, 5)
              sage: i^2
              -3
              sage: j^2
              5
          """
          return self._a, self._b
  
      def __eq__(self, other):
          """
          Compare self and other.
  
          EXAMPLES::
  
              sage: QuaternionAlgebra(-1,-7) == QuaternionAlgebra(-1,-7)
              True
              sage: QuaternionAlgebra(-1,-7) == QuaternionAlgebra(-1,-5)
              False
          """
          if not isinstance(other, QuaternionAlgebra_abstract):
              return False
          return (self.base_ring() == other.base_ring() and
                  (self._a, self._b) == (other._a, other._b))
  
      def __ne__(self, other):
          """
          Compare self and other.
  
          EXAMPLES::
  
              sage: QuaternionAlgebra(-1,-7) != QuaternionAlgebra(-1,-7)
              False
              sage: QuaternionAlgebra(-1,-7) != QuaternionAlgebra(-1,-5)
              True
          """
          return not self.__eq__(other)
  
      def __hash__(self):
          """
          Compute the hash of ``self``.
  
          EXAMPLES::
  
              sage: h1 = hash(QuaternionAlgebra(-1,-7))
              sage: h2 = hash(QuaternionAlgebra(-1,-7))
              sage: h3 = hash(QuaternionAlgebra(-1,-5))
              sage: h1  == h2 and h1 != h3
              True
          """
          return hash((self.base_ring(), self._a, self._b))
  
      def gen(self, i=0):
          """
          Return the `i^{th}` generator of ``self``.
  
          INPUT:
  
          - ``i`` - integer (optional, default 0)
  
          EXAMPLES::
  
              sage: Q.<ii,jj,kk> = QuaternionAlgebra(QQ,-1,-2); Q
              Quaternion Algebra (-1, -2) with base ring Rational Field
              sage: Q.gen(0)
              ii
              sage: Q.gen(1)
              jj
              sage: Q.gen(2)
              kk
              sage: Q.gens()
              [ii, jj, kk]
          """
          return self._gens[i]
  
      def _repr_(self):
          """
          Print representation.
  
          TESTS::
  
              sage: Q.<i,j,k> = QuaternionAlgebra(QQ,-5,-2)
              sage: type(Q)
              <class 'sage.algebras.quatalg.quaternion_algebra.QuaternionAlgebra_ab_with_category'>
              sage: Q._repr_()
              'Quaternion Algebra (-5, -2) with base ring Rational Field'
              sage: Q
              Quaternion Algebra (-5, -2) with base ring Rational Field
              sage: print(Q)
              Quaternion Algebra (-5, -2) with base ring Rational Field
              sage: str(Q)
              'Quaternion Algebra (-5, -2) with base ring Rational Field'
          """
          return "Quaternion Algebra (%r, %r) with base ring %s" % (self._a, self._b, self.base_ring())
  
      def inner_product_matrix(self):
          """
          Return the inner product matrix associated to ``self``, i.e. the
          Gram matrix of the reduced norm as a quadratic form on ``self``.
          The standard basis `1`, `i`, `j`, `k` is orthogonal, so this matrix
          is just the diagonal matrix with diagonal entries `1`, `a`, `b`, `ab`.
  
          EXAMPLES::
  
              sage: Q.<i,j,k> = QuaternionAlgebra(-5,-19)
              sage: Q.inner_product_matrix()
              [  2   0   0   0]
              [  0  10   0   0]
              [  0   0  38   0]
              [  0   0   0 190]
  
              sage: R.<a,b> = QQ[]; Q.<i,j,k> = QuaternionAlgebra(Frac(R),a,b)
              sage: Q.inner_product_matrix()
              [    2     0     0     0]
              [    0  -2*a     0     0]
              [    0     0  -2*b     0]
              [    0     0     0 2*a*b]
          """
          a, b = self._a, self._b
          return diagonal_matrix(self.base_ring(), [2, -2*a, -2*b, 2*a*b])
  
      @cached_method
      def discriminant(self):
          """
          Given a quaternion algebra `A` defined over a number field,
          return the discriminant of `A`, i.e. the
          product of the ramified primes of `A`.
  
          EXAMPLES::
  
              sage: QuaternionAlgebra(210,-22).discriminant()
              210
              sage: QuaternionAlgebra(19).discriminant()
              19
  
              sage: F.<a> = NumberField(x^2-x-1)
              sage: B.<i,j,k> = QuaternionAlgebra(F, 2*a,F(-1))
              sage: B.discriminant()
              Fractional ideal (2)
  
              sage: QuaternionAlgebra(QQ[sqrt(2)],3,19).discriminant()
              Fractional ideal (1)
          """
          if not is_RationalField(self.base_ring()):
              try:
                  F = self.base_ring()
                  return F.hilbert_conductor(self._a, self._b)
              except NotImplementedError:
                  raise ValueError("base field must be rational numbers or number field")
          else:
              return hilbert_conductor(self._a, self._b)
  
      def ramified_primes(self):
          """
          Return the primes that ramify in this quaternion algebra.
  
          Currently only implemented over the rational numbers.
  
          EXAMPLES::
  
              sage: QuaternionAlgebra(QQ, -1, -1).ramified_primes()
              [2]
          """
          # TODO: more examples
          return [f[0] for f in factor(self.discriminant())]
  
      def _magma_init_(self, magma):
          """
          Return Magma version of this quaternion algebra.
  
          EXAMPLES::
  
              sage: Q = QuaternionAlgebra(-1,-1); Q
              Quaternion Algebra (-1, -1) with base ring Rational Field
              sage: Q._magma_init_(magma)                                  # optional - magma
              'QuaternionAlgebra(_sage_[...],-1/1,-1/1)'
              sage: A = magma(Q); A                                        # optional - magma
              Quaternion Algebra with base ring Rational Field, defined by i^2 = -1, j^2 = -1
              sage: A.RamifiedPlaces()                                     # optional - magma
              [
              Ideal of Integer Ring generated by 2
              ]
  
          A more complicated example involving a quaternion algebra over a number field::
  
              sage: K.<a> = QQ[sqrt(2)]; Q = QuaternionAlgebra(K,-1,a); Q
              Quaternion Algebra (-1, sqrt2) with base ring Number Field in sqrt2 with defining polynomial x^2 - 2 with sqrt2 = 1.414213562373095?
              sage: magma(Q)                                              # optional - magma
              Quaternion Algebra with base ring Number Field with defining polynomial x^2 - 2 over the Rational Field, defined by i^2 = -1, j^2 = sqrt2
              sage: Q._magma_init_(magma)                                 # optional - magma
              'QuaternionAlgebra(_sage_[...],(_sage_[...]![-1, 0]),(_sage_[...]![0, 1]))'
          """
          R = magma(self.base_ring())
          return 'QuaternionAlgebra(%s,%s,%s)' % (R.name(),
                                                  self._a._magma_init_(magma),
                                                  self._b._magma_init_(magma))
  
      def quaternion_order(self, basis, check=True):
          """
          Return the order of this quaternion order with given basis.
  
          INPUT:
  
          - ``basis`` - list of 4 elements of ``self``
          - ``check`` - bool (default: ``True``)
  
          EXAMPLES::
  
              sage: Q.<i,j,k> = QuaternionAlgebra(-11,-1)
              sage: Q.quaternion_order([1,i,j,k])
              Order of Quaternion Algebra (-11, -1) with base ring Rational Field with basis (1, i, j, k)
  
          We test out ``check=False``::
  
              sage: Q.quaternion_order([1,i,j,k], check=False)
              Order of Quaternion Algebra (-11, -1) with base ring Rational Field with basis (1, i, j, k)
              sage: Q.quaternion_order([i,j,k], check=False)
              Order of Quaternion Algebra (-11, -1) with base ring Rational Field with basis (i, j, k)
          """
          return QuaternionOrder(self, basis, check=check)
  
      def ideal(self, gens, left_order=None, right_order=None, check=True, **kwds):
          r"""
          Return the quaternion ideal with given gens over `\ZZ`.
  
          Neither a left or right order structure need be specified.
  
          INPUT:
  
          - ``gens`` -- a list of elements of this quaternion order
  
          - ``check`` -- bool (default: ``True``)
  
          - ``left_order`` -- a quaternion order or ``None``
  
          - ``right_order`` -- a quaternion order or ``None``
  
          EXAMPLES::
  
              sage: R = QuaternionAlgebra(-11,-1)
              sage: R.ideal([2*a for a in R.basis()])
              Fractional ideal (2, 2*i, 2*j, 2*k)
          """
          gens = [self(g) for g in gens]  # coerce integers etc. into quaternions
          if self.base_ring() == QQ:
              return QuaternionFractionalIdeal_rational(self, gens, left_order=left_order, right_order=right_order, check=check)
          else:
              raise NotImplementedError("ideal only implemented for quaternion algebras over QQ")
  
      @cached_method
      def modp_splitting_data(self, p):
          r"""
          Return mod `p` splitting data for this quaternion algebra at
          the unramified prime `p`.
  
          This is `2\times 2`
          matrices `I`, `J`, `K` over the finite field `\GF{p}` such that if
          the quaternion algebra has generators `i, j, k`, then `I^2 =
          i^2`, `J^2 = j^2`, `IJ=K` and `IJ=-JI`.
  
          .. NOTE::
  
             Currently only implemented when `p` is odd and the base
             ring is `\QQ`.
  
          INPUT:
  
          - `p` -- unramified odd prime
  
          OUTPUT:
  
          - 2-tuple of matrices over finite field
  
          EXAMPLES::
  
              sage: Q = QuaternionAlgebra(-15, -19)
              sage: Q.modp_splitting_data(7)
              (
              [0 6]  [6 1]  [6 6]
              [1 0], [1 1], [6 1]
              )
              sage: Q.modp_splitting_data(next_prime(10^5))
              (
              [    0 99988]  [97311     4]  [99999 59623]
              [    1     0], [13334  2692], [97311     4]
              )
              sage: I,J,K = Q.modp_splitting_data(23)
              sage: I
              [0 8]
              [1 0]
              sage: I^2
              [8 0]
              [0 8]
              sage: J
              [19  2]
              [17  4]
              sage: J^2
              [4 0]
              [0 4]
              sage: I*J == -J*I
              True
              sage: I*J == K
              True
  
          The following is a good test because of the asserts in the code::
  
              sage: v = [Q.modp_splitting_data(p) for p in primes(20,1000)]
  
          Proper error handling::
  
              sage: Q.modp_splitting_data(5)
              Traceback (most recent call last):
              ...
              NotImplementedError: algorithm for computing local splittings not implemented in general (currently require the first invariant to be coprime to p)
  
              sage: Q.modp_splitting_data(2)
              Traceback (most recent call last):
              ...
              NotImplementedError: p must be odd
          """
          if self.base_ring() != QQ:
              raise NotImplementedError("must be rational quaternion algebra")
          p = ZZ(p)
          if not p.is_prime():
              raise ValueError("p (=%s) must be prime" % p)
          if p == 2:
              raise NotImplementedError("p must be odd")
          if self.discriminant() % p == 0:
              raise ValueError("p (=%s) must be an unramified prime" % p)
  
          i, j, _ = self.gens()
          F = GF(p)
          i2 = F(i * i)
          j2 = F(j * j)
  
          M = MatrixSpace(F, 2)
          I = M([0, i2, 1, 0])
          if i2 == 0:
              raise NotImplementedError("algorithm for computing local splittings not implemented in general (currently require the first invariant to be coprime to p)")
          i2inv = ~i2
          a = None
          for b in list(F):
              if not b:
                  continue
              c = j2 + i2inv * b*b
              if c.is_square():
                  a = -c.sqrt()
                  break
  
          if a is None:
              # do a fallback search, maybe needed in char 3 sometimes.
              for J in M:
                  K = I*J
                  if J*J == j2 and K == -J*I:
                      return I, J, K
  
          J = M([a, b, (j2 - a * a) / b, -a])
          K = I * J
          assert K == -J * I, "bug in that I,J don't skew commute"
          return I, J, K
  
      def modp_splitting_map(self, p):
          r"""
          Return Python map from the (`p`-integral) quaternion algebra to
          the set of `2\times 2` matrices over `\GF{p}`.
  
          INPUT:
  
          - `p` -- prime number
  
          EXAMPLES::
  
              sage: Q.<i,j,k> = QuaternionAlgebra(-1, -7)
              sage: f = Q.modp_splitting_map(13)
              sage: a = 2+i-j+3*k; b = 7+2*i-4*j+k
              sage: f(a*b)
              [12  3]
              [10  5]
              sage: f(a)*f(b)
              [12  3]
              [10  5]
          """
          I, J, K = self.modp_splitting_data(p)
          F = I.base_ring()
  
          def phi(q):
              v = [F(a) for a in q.coefficient_tuple()]
              return v[0] + I*v[1] + J*v[2] + K*v[3]
          return phi
  
  
  ############################################################
  # Unpickling
  ############################################################
  def unpickle_QuaternionAlgebra_v0(*key):
      """
      The 0th version of pickling for quaternion algebras.
  
      EXAMPLES::
  
          sage: Q = QuaternionAlgebra(-5,-19)
          sage: t = (QQ, -5, -19, ('i', 'j', 'k'))
          sage: sage.algebras.quatalg.quaternion_algebra.unpickle_QuaternionAlgebra_v0(*t)
          Quaternion Algebra (-5, -19) with base ring Rational Field
          sage: loads(dumps(Q)) == Q
          True
          sage: loads(dumps(Q)) is Q
          True
      """
      return QuaternionAlgebra(*key)
  
  
  class QuaternionOrder(Parent):
      """
      An order in a quaternion algebra.
  
      EXAMPLES::
  
          sage: QuaternionAlgebra(-1,-7).maximal_order()
          Order of Quaternion Algebra (-1, -7) with base ring Rational Field with basis (1/2 + 1/2*j, 1/2*i + 1/2*k, j, k)
          sage: type(QuaternionAlgebra(-1,-7).maximal_order())
          <class 'sage.algebras.quatalg.quaternion_algebra.QuaternionOrder_with_category'>
      """
      def __init__(self, A, basis, check=True):
          """
          INPUT:
  
          - ``A`` - a quaternion algebra
          - ``basis`` - list of 4 integral quaternions in ``A``
          - ``check`` - whether to do type and other consistency checks
  
          .. WARNING::
  
              Currently most methods silently assume that the ``A.base_ring()``
              is ``QQ``.
  
          EXAMPLES::
  
              sage: A.<i,j,k> = QuaternionAlgebra(-3,-5)
              sage: sage.algebras.quatalg.quaternion_algebra.QuaternionOrder(A, [1,i,j,k])
              Order of Quaternion Algebra (-3, -5) with base ring Rational Field with basis (1, i, j, k)
              sage: R = sage.algebras.quatalg.quaternion_algebra.QuaternionOrder(A, [1,2*i,2*j,2*k]); R
              Order of Quaternion Algebra (-3, -5) with base ring Rational Field with basis (1, 2*i, 2*j, 2*k)
              sage: type(R)
              <class 'sage.algebras.quatalg.quaternion_algebra.QuaternionOrder_with_category'>
  
              Over QQ and number fields it is checked whether the given
              basis actually gives an order (as a module over the maximal order):
  
              sage: A.<i,j,k> = QuaternionAlgebra(-1,-1)
              sage: A.quaternion_order([1,i,j,i-j])
              Traceback (most recent call last):
              ...
              ValueError: basis must have rank 4
              sage: A.quaternion_order([2,i,j,k])
              Traceback (most recent call last):
              ...
              ValueError: lattice must contain 1
              sage: A.quaternion_order([1,i/2,j/2,k/2])
              Traceback (most recent call last):
              ...
              ValueError: given lattice must be a ring
  
              sage: K = QuadraticField(10)
              sage: A.<i,j,k> = QuaternionAlgebra(K,-1,-1)
              sage: A.quaternion_order([1,i,j,k])
              Order of Quaternion Algebra (-1, -1) with base ring Number Field in a with defining polynomial x^2 - 10 with a = 3.162277660168380? with basis (1, i, j, k)
              sage: A.quaternion_order([1,i/2,j,k])
              Traceback (most recent call last):
              ...
              ValueError: given lattice must be a ring
  
          TESTS::
  
              sage: TestSuite(R).run()
          """
          if check:
              # right data type
              if not isinstance(basis, (list, tuple)):
                  raise TypeError("basis must be a list or tuple")
              # right length
              if len(basis) != 4:
                  raise ValueError("basis must have length 4")
              # coerce to common parent
              basis = tuple([A(x) for x in basis])
  
              # has rank 4
              V = A.base_ring()**4
              if V.span([V(x.coefficient_tuple()) for x in basis]).dimension() != 4:
                  raise ValueError("basis must have rank 4")
  
              # The additional checks will work over QQ and over number fields,
              # but we can't actually do much with an order defined over a number
              # field
  
              if A.base_ring() == QQ:     # fast code over QQ
                  M = matrix(QQ, 4, 4, [x.coefficient_tuple() for x in basis])
                  v = M.solve_left(V([1, 0, 0, 0]))
  
                  if v.denominator() != 1:
                      raise ValueError("lattice must contain 1")
  
                  # check if multiplicatively closed
                  M1 = basis_for_quaternion_lattice(basis)
                  M2 = basis_for_quaternion_lattice(list(basis) + [x * y for x in basis for y in basis])
                  if M1 != M2:
                      raise ValueError("given lattice must be a ring")
  
              if A.base_ring() != QQ:     # slow code over number fields (should eventually use PARI's nfhnf)
                  O = None
                  try:
                      O = A.base_ring().maximal_order()
                  except AttributeError:
                      pass
  
                  if O:
                      M = matrix(A.base_ring(), 4, 4, [x.coefficient_tuple()
                                                       for x in basis])
                      v = M.solve_left(V([1, 0, 0, 0]))
  
                      if any(a not in O for a in v):
                          raise ValueError("lattice must contain 1")
  
                      # check if multiplicatively closed
                      Y = matrix(QQ, 16, 4, [(x*y).coefficient_tuple()
                                             for x in basis for y in basis])
                      X = M.solve_left(Y)
                      if any(a not in O for x in X for a in x):
                          raise ValueError("given lattice must be a ring")
  
          self.__basis = tuple(basis)
          self.__quaternion_algebra = A
          Parent.__init__(self, base=ZZ, facade=(A,),
                          category=Algebras(ZZ).Facade().FiniteDimensional())
  
      def _element_constructor_(self, x):
          """
          Construct an element of this quaternion order from ``x``,
          or throw an error if ``x`` is not contained in the order.
  
          EXAMPLES::
  
              sage: Q.<i,j,k> = QuaternionAlgebra(-1,-19)
              sage: O = Q.quaternion_order([1,i,j,k])
              sage: O(1+i)
              1 + i
              sage: O(1/2)
              Traceback (most recent call last):
              ...
              TypeError: 1/2 does not lie in Order of Quaternion Algebra (-1, -19)
              with base ring Rational Field with basis (1, i, j, k)
  
          TESTS:
  
          Test for :trac:`32364`::
  
              sage: 1/5 in O
              False
              sage: j/2 in O
              False
  
          """
          y = self.quaternion_algebra()(x)
          if y not in self.unit_ideal():
              raise TypeError(f'{x!r} does not lie in {self!r}')
          return y
  
      def one(self):
          """
          Return the multiplicative unit of this quaternion order.
  
          EXAMPLES::
  
              sage: QuaternionAlgebra(-1,-7).maximal_order().one()
              1
          """
          return self.quaternion_algebra().one()
  
      def gens(self):
          """
          Return generators for self.
  
          EXAMPLES::
  
              sage: QuaternionAlgebra(-1,-7).maximal_order().gens()
              (1/2 + 1/2*j, 1/2*i + 1/2*k, j, k)
          """
          return self.__basis
  
      def ngens(self):
          """
          Return the number of generators (which is 4).
  
          EXAMPLES::
  
              sage: QuaternionAlgebra(-1,-7).maximal_order().ngens()
              4
          """
          return 4
  
      def gen(self, n):
          """
          Return the n-th generator.
  
          INPUT:
  
          - ``n`` - an integer between 0 and 3, inclusive.
  
          EXAMPLES::
  
              sage: R = QuaternionAlgebra(-11,-1).maximal_order(); R
              Order of Quaternion Algebra (-11, -1) with base ring Rational Field with basis (1/2 + 1/2*i, 1/2*j - 1/2*k, i, -k)
              sage: R.gen(0)
              1/2 + 1/2*i
              sage: R.gen(1)
              1/2*j - 1/2*k
              sage: R.gen(2)
              i
              sage: R.gen(3)
              -k
          """
          return self.__basis[n]
  
      def __eq__(self, R):
          """
          Compare orders self and other.
  
          EXAMPLES::
  
              sage: R = QuaternionAlgebra(-11,-1).maximal_order()
              sage: R == R                       # indirect doctest
              True
              sage: R == QuaternionAlgebra(-1,-1).maximal_order()
              False
              sage: R == 5
              False
              sage: Q.<i,j,k> = QuaternionAlgebra(-1,-19)
  
          Orders can be equal even if they are defined by different
          bases (see :trac:`32245`)::
  
              sage: Q.quaternion_order([1,-i,k,j+i*7]) == Q.quaternion_order([1,i,j,k])
              True
          """
          if not isinstance(R, QuaternionOrder):
              return False
          return (self.__quaternion_algebra == R.__quaternion_algebra and
                  self.unit_ideal() == R.unit_ideal())
  
      def __ne__(self, other):
          """
          Compare orders self and other.
  
          Two orders are equal if they
          have the same basis and are in the same quaternion algebra.
  
          EXAMPLES::
  
              sage: R = QuaternionAlgebra(-11,-1).maximal_order()
              sage: R != R                       # indirect doctest
              False
              sage: R != QuaternionAlgebra(-1,-1).maximal_order()
              True
          """
          return not self.__eq__(other)
  
      def __hash__(self):
          """
          Compute the hash of ``self``.
  
          EXAMPLES::
  
              sage: h1 = hash(QuaternionAlgebra(-1,-7).maximal_order())
              sage: h2 = hash(QuaternionAlgebra(-1,-7).maximal_order())
              sage: h3 = hash(QuaternionAlgebra(-1,-5).maximal_order())
              sage: h1  == h2 and h1 != h3
              True
          """
          return hash((self.__quaternion_algebra, self.__basis))
  
      def basis(self):
          """
          Return fix choice of basis for this quaternion order.
  
          EXAMPLES::
  
              sage: QuaternionAlgebra(-11,-1).maximal_order().basis()
              (1/2 + 1/2*i, 1/2*j - 1/2*k, i, -k)
          """
          return self.__basis
  
      def quaternion_algebra(self):
          """
          Return ambient quaternion algebra that contains this quaternion order.
  
          EXAMPLES::
  
              sage: QuaternionAlgebra(-11,-1).maximal_order().quaternion_algebra()
              Quaternion Algebra (-11, -1) with base ring Rational Field
          """
          return self.__quaternion_algebra
  
      def _repr_(self):
          """
          Return string representation of this order.
  
          EXAMPLES::
  
              sage: QuaternionAlgebra(-11,-1).maximal_order()._repr_()
              'Order of Quaternion Algebra (-11, -1) with base ring Rational Field with basis (1/2 + 1/2*i, 1/2*j - 1/2*k, i, -k)'
              sage: QuaternionAlgebra(-11,-1).maximal_order()
              Order of Quaternion Algebra (-11, -1) with base ring Rational Field with basis (1/2 + 1/2*i, 1/2*j - 1/2*k, i, -k)
          """
          return 'Order of %s with basis %s' % (self.quaternion_algebra(), self.basis())
  
      def random_element(self, *args, **kwds):
          """
          Return a random element of this order.
  
          The args and kwds are passed to the random_element method of
          the integer ring, and we return an element of the form
  
          .. MATH::
  
              ae_1 + be_2 + ce_3 + de_4
  
          where `e_1`, ..., `e_4` are the basis of this order and `a`,
          `b`, `c`, `d` are random integers.
  
          EXAMPLES::
  
              sage: QuaternionAlgebra(-11,-1).maximal_order().random_element()  # random
              -4 - 4*i + j - k
              sage: QuaternionAlgebra(-11,-1).maximal_order().random_element(-10,10)  # random
              -9/2 - 7/2*i - 7/2*j - 3/2*k
          """
          return sum(ZZ.random_element(*args, **kwds) * b for b in self.basis())
  
      def intersection(self, other):
          """
          Return the intersection of this order with other.
  
          INPUT:
  
          - ``other`` - a quaternion order in the same ambient quaternion algebra
  
          OUTPUT: a quaternion order
  
          EXAMPLES::
  
              sage: R = QuaternionAlgebra(-11,-1).maximal_order()
              sage: R.intersection(R)
              Order of Quaternion Algebra (-11, -1) with base ring Rational Field with basis (1/2 + 1/2*i, i, 1/2*j + 1/2*k, k)
  
          We intersect various orders in the quaternion algebra ramified at 11::
  
              sage: B = BrandtModule(11,3)
              sage: R = B.maximal_order(); S = B.order_of_level_N()
              sage: R.intersection(S)
              Order of Quaternion Algebra (-1, -11) with base ring Rational Field with basis (1/2 + 1/2*j, 1/2*i + 5/2*k, j, 3*k)
              sage: R.intersection(S) == S
              True
              sage: B = BrandtModule(11,5)
              sage: T = B.order_of_level_N()
              sage: S.intersection(T)
              Order of Quaternion Algebra (-1, -11) with base ring Rational Field with basis (1/2 + 1/2*j, 1/2*i + 23/2*k, j, 15*k)
          """
          if not isinstance(other, QuaternionOrder):
              raise TypeError("other must be a QuaternionOrder")
  
          A = self.quaternion_algebra()
          if other.quaternion_algebra() != A:
              raise ValueError("self and other must be in the same ambient quaternion algebra")
  
          V = A.base_ring()**4
  
          B = V.span([V(list(g)) for g in self.basis()], ZZ)
          C = V.span([V(list(g)) for g in other.basis()], ZZ)
  
          # todo -- A(list(e)) could be A(e)
          return QuaternionOrder(A, [A(list(e)) for e in B.intersection(C).basis()])
  
      @cached_method
      def free_module(self):
          r"""
          Return the free `\ZZ`-module that corresponds to this order
          inside the vector space corresponding to the ambient
          quaternion algebra.
  
          OUTPUT:
  
          A free `\ZZ`-module of rank 4.
  
          EXAMPLES::
  
              sage: R = QuaternionAlgebra(-11,-1).maximal_order()
              sage: R.basis()
              (1/2 + 1/2*i, 1/2*j - 1/2*k, i, -k)
              sage: R.free_module()
              Free module of degree 4 and rank 4 over Integer Ring
              Echelon basis matrix:
              [1/2 1/2   0   0]
              [  0   1   0   0]
              [  0   0 1/2 1/2]
              [  0   0   0   1]
          """
          V = self.quaternion_algebra().base_ring()**4
          return V.span([V(list(g)) for g in self.basis()], ZZ)
  
      def discriminant(self):
          r"""
          Return the discriminant of this order.
  
          This is defined as
          `\sqrt{ det ( Tr(e_i \bar{e}_j ) ) }`, where `\{e_i\}` is the
          basis of the order.
  
          OUTPUT: rational number
  
          EXAMPLES::
  
              sage: QuaternionAlgebra(-11,-1).maximal_order().discriminant()
              11
              sage: S = BrandtModule(11,5).order_of_level_N()
              sage: S.discriminant()
              55
              sage: type(S.discriminant())
              <... 'sage.rings.rational.Rational'>
          """
          L = []
          for d in self.basis():
              MM = []
              for e in self.basis():
                  MM.append(d.pair(e))
              L.append(MM)
  
          return (MatrixSpace(QQ, 4, 4)(L)).determinant().sqrt()
  
      def left_ideal(self, gens, check=True):
          r"""
          Return the ideal with given gens over `\ZZ`.
  
          INPUT:
  
          - ``gens`` -- a list of elements of this quaternion order
  
          - ``check`` -- bool (default: ``True``)
  
          EXAMPLES::
  
              sage: R = QuaternionAlgebra(-11,-1).maximal_order()
              sage: R.left_ideal([2*a for a in R.basis()])
              Fractional ideal (1 + i, 2*i, j + k, 2*k)
          """
          if self.base_ring() == ZZ:
              return QuaternionFractionalIdeal_rational(self.quaternion_algebra(), gens, left_order=self, check=check)
          else:
              raise NotImplementedError("ideal only implemented for quaternion algebras over QQ")
  
      def right_ideal(self, gens, check=True):
          r"""
          Return the ideal with given gens over `\ZZ`.
  
          INPUT:
  
          - ``gens`` -- a list of elements of this quaternion order
  
          - ``check`` -- bool (default: ``True``)
  
          EXAMPLES::
  
              sage: R = QuaternionAlgebra(-11,-1).maximal_order()
              sage: R.right_ideal([2*a for a in R.basis()])
              Fractional ideal (1 + i, 2*i, j + k, 2*k)
          """
          if self.base_ring() == ZZ:
              return QuaternionFractionalIdeal_rational(self.quaternion_algebra(), gens, right_order=self, check=check)
          else:
              raise NotImplementedError("ideal only implemented for quaternion algebras over QQ")
  
      @cached_method
      def unit_ideal(self):
          """
          Return the unit ideal in this quaternion order.
  
          EXAMPLES::
  
              sage: R = QuaternionAlgebra(-11,-1).maximal_order()
              sage: I = R.unit_ideal(); I
              Fractional ideal (1/2 + 1/2*i, 1/2*j - 1/2*k, i, -k)
          """
          if self.base_ring() == ZZ:
              return QuaternionFractionalIdeal_rational(self.quaternion_algebra(), self.basis(), left_order=self, right_order=self, check=False)
          else:
              raise NotImplementedError("ideal only implemented for quaternion algebras over QQ")
  
      def __mul__(self, other):
          """
          Every order equals its own unit ideal. Overload ideal multiplication
          and scaling to orders.
  
          EXAMPLES::
  
              sage: Q.<i,j,k> = QuaternionAlgebra(-1,-11)
              sage: O = Q.maximal_order()
              sage: I = O*j; I
              Fractional ideal (-11/2 + 1/2*j, -11/2*i + 1/2*k, -11, -11*i)
          """
          return self.unit_ideal() * other
  
      def __rmul__(self, other):
          return other * self.unit_ideal()
  
      def __add__(self, other):
          """
          Every order equals its own unit ideal. Overload ideal addition
          to orders.
  
          EXAMPLES::
  
              sage: Q.<i,j,k> = QuaternionAlgebra(-1,-11)
              sage: O = Q.maximal_order()
              sage: I = O + O*((j-3)/5); I
              Fractional ideal (1/10 + 3/10*j, 1/10*i + 3/10*k, j, k)
          """
          return self.unit_ideal() + other
  
      def quadratic_form(self):
          """
          Return the normalized quadratic form associated to this quaternion order.
  
          OUTPUT: quadratic form
  
          EXAMPLES::
  
              sage: R = BrandtModule(11,13).order_of_level_N()
              sage: Q = R.quadratic_form(); Q
              Quadratic form in 4 variables over Rational Field with coefficients:
              [ 14 253 55 286 ]
              [ * 1455 506 3289 ]
              [ * * 55 572 ]
              [ * * * 1859 ]
              sage: Q.theta_series(10)
              1 + 2*q + 2*q^4 + 4*q^6 + 4*q^8 + 2*q^9 + O(q^10)
          """
          return self.unit_ideal().quadratic_form()
  
      def ternary_quadratic_form(self, include_basis=False):
          r"""
          Return the ternary quadratic form associated to this order.
  
          INPUT:
  
          - ``include_basis`` -- bool (default: False), if True also
            return a basis for the dimension 3 subspace `G`
  
          OUTPUT:
  
          - QuadraticForm
  
          - optional basis for dimension 3 subspace
  
          This function computes the positive definition quadratic form
          obtained by letting G be the trace zero subspace of `\ZZ` +
          2* ``self``, which has rank 3, and restricting the pairing
          :meth:`QuaternionAlgebraElement_abstract.pair`::
  
             (x,y) = (x.conjugate()*y).reduced_trace()
  
          to `G`.
  
          APPLICATIONS: Ternary quadratic forms associated to an order
          in a rational quaternion algebra are useful in computing with
          Gross points, in decided whether quaternion orders have
          embeddings from orders in quadratic imaginary fields, and in
          computing elements of the Kohnen plus subspace of modular
          forms of weight 3/2.
  
          EXAMPLES::
  
              sage: R = BrandtModule(11,13).order_of_level_N()
              sage: Q = R.ternary_quadratic_form(); Q
              Quadratic form in 3 variables over Rational Field with coefficients:
              [ 5820 1012 13156 ]
              [ * 55 1144 ]
              [ * * 7436 ]
              sage: factor(Q.disc())
              2^4 * 11^2 * 13^2
  
          The following theta series is a modular form of weight 3/2 and level 4*11*13::
  
              sage: Q.theta_series(100)
              1 + 2*q^23 + 2*q^55 + 2*q^56 + 2*q^75 + 4*q^92 + O(q^100)
          """
          if self.base_ring() != ZZ:
              raise NotImplementedError("ternary quadratic form of order only implemented for quaternion algebras over QQ")
  
          Q = self.quaternion_algebra()
          # 2*R + ZZ
          twoR = self.free_module().scale(2)
          Z = twoR.span([Q(1).coefficient_tuple()], ZZ)
          S = twoR + Z
          # Now we intersect with the trace 0 submodule
          v = [b.reduced_trace() for b in Q.basis()]
          M = matrix(QQ, 4, 1, v)
          tr0 = M.kernel()
          G = tr0.intersection(S)
          B = [Q(a) for a in G.basis()]
          m = matrix(QQ, [[x.pair(y) for x in B] for y in B])
          from sage.quadratic_forms.quadratic_form import QuadraticForm
          Q = QuadraticForm(m)
          if include_basis:
              return Q, B
          else:
              return Q
  
  
  class QuaternionFractionalIdeal(Ideal_fractional):
      pass
  
  
  class QuaternionFractionalIdeal_rational(QuaternionFractionalIdeal):
      """
      A fractional ideal in a rational quaternion algebra.
  
      INPUT:
  
      - ``left_order`` -- a quaternion order or ``None``
  
      - ``right_order`` -- a quaternion order or ``None``
  
      - ``basis`` -- tuple of length 4 of elements in of ambient
        quaternion algebra whose `\\ZZ`-span is an ideal
  
      - ``check`` -- bool (default: ``True``); if ``False``, do no type
        checking.
      """
      def __init__(self, Q, basis, left_order=None, right_order=None, check=True):
          """
          EXAMPLES::
  
              sage: R = QuaternionAlgebra(-11,-1).maximal_order()
              sage: R.right_ideal(R.basis())
              Fractional ideal (1/2 + 1/2*i, i, 1/2*j + 1/2*k, k)
              sage: R.right_ideal(tuple(R.basis()), check=False)
              Fractional ideal (1/2 + 1/2*i, 1/2*j - 1/2*k, i, -k)
  
          TESTS::
  
              sage: QuaternionAlgebra(-11,-1).maximal_order().unit_ideal().gens()
              (1/2 + 1/2*i, 1/2*j - 1/2*k, i, -k)
  
          Check that :trac:`31582` is fixed::
  
              sage: BrandtModule(23).right_ideals()[0].parent()
              Monoid of ideals of Quaternion Algebra (-1, -23) with base ring Rational Field
          """
          if check:
              if left_order is not None and not isinstance(left_order, QuaternionOrder):
                  raise TypeError("left_order must be a quaternion order or None")
              if right_order is not None and not isinstance(right_order, QuaternionOrder):
                  raise TypeError("right_order must be a quaternion order or None")
              if not isinstance(basis, (list, tuple)):
                  raise TypeError("basis must be a list or tuple")
              basis = tuple([Q(v) for v in
                             (QQ**4).span([Q(v).coefficient_tuple() for v in basis], ZZ).basis()])
          self.__left_order = left_order
          self.__right_order = right_order
          Ideal_fractional.__init__(self, Q, basis)
  
      def scale(self, alpha, left=False):
          r"""
          Scale the fractional ideal ``self`` by multiplying the basis
          by ``alpha``.
  
          INPUT:
  
          - `\alpha` -- element of quaternion algebra
  
          - ``left`` -- bool (default: False); if true multiply
            `\alpha` on the left, otherwise multiply `\alpha` on the right
  
          OUTPUT:
  
          - a new fractional ideal
  
          EXAMPLES::
  
              sage: B = BrandtModule(5,37); I = B.right_ideals()[0]; i,j,k = B.quaternion_algebra().gens(); I
              Fractional ideal (2 + 2*j + 106*k, i + 2*j + 105*k, 4*j + 64*k, 148*k)
              sage: I.scale(i)
              Fractional ideal (2*i + 212*j - 2*k, -2 + 210*j - 2*k, 128*j - 4*k, 296*j)
              sage: I.scale(i, left=True)
              Fractional ideal (2*i - 212*j + 2*k, -2 - 210*j + 2*k, -128*j + 4*k, -296*j)
              sage: I.scale(i, left=False)
              Fractional ideal (2*i + 212*j - 2*k, -2 + 210*j - 2*k, 128*j - 4*k, 296*j)
              sage: i * I.gens()[0]
              2*i - 212*j + 2*k
              sage: I.gens()[0] * i
              2*i + 212*j - 2*k
  
          TESTS:
  
          Scaling by `1` should not change anything (see :trac:`32245`)::
  
              sage: I.scale(1) == I
              True
  
          Check that :trac:`32726` is fixed::
  
              sage: Q.<i,j,k> = QuaternionAlgebra(-1,-19)
              sage: I = Q.ideal([1,i,j,k])
              sage: _ = I.left_order(), I.right_order()   # cache
              sage: span = lambda L: L.basis_matrix().row_module(ZZ)
              sage: for left in (True,False):
              ....:     J = I.scale(1+i+j+k, left=left)
              ....:     Ol, Or = J.left_order(), J.right_order()
              ....:     [
              ....:       span(Ol.unit_ideal() * J) <= span(J),
              ....:       span(J * Or.unit_ideal()) <= span(J),
              ....:     ]
              [True, True]
              [True, True]
          """
          Q = self.quaternion_algebra()
          alpha = Q(alpha)
          if left:
              gens = [alpha * b for b in self.basis()]
          else:
              gens = [b * alpha for b in self.basis()]
          left_order = self.__left_order if alpha in QQ or not left else None
          right_order = self.__right_order if alpha in QQ or left else None
          return Q.ideal(gens, check=False,
                         left_order=left_order, right_order=right_order)
  
      def quaternion_algebra(self):
          """
          Return the ambient quaternion algebra that contains this fractional ideal.
  
          OUTPUT: a quaternion algebra
  
          EXAMPLES::
  
              sage: I = BrandtModule(3,5).right_ideals()[1]; I
              Fractional ideal (2 + 6*j + 4*k, 2*i + 4*j + 34*k, 8*j + 32*k, 40*k)
              sage: I.quaternion_algebra()
              Quaternion Algebra (-1, -3) with base ring Rational Field
          """
          # TODO: when the ring() method is removed from this class, the
          # following line can be changed to self.ring().
          return Ideal_fractional.ring(self)
  
      def _compute_order(self, side='left'):
          r"""
          Used internally to compute either the left or right order
          associated to an ideal in a quaternion algebra.  If
          action='right', compute the left order, and if action='left'
          compute the right order.
  
          INPUT:
  
          - ``side`` -- 'left' or 'right'
  
          EXAMPLES::
  
              sage: R.<i,j,k> = QuaternionAlgebra(-1,-11)
              sage: I = R.ideal([2 + 2*j + 140*k, 2*i + 4*j + 150*k, 8*j + 104*k, 152*k])
              sage: Ol = I._compute_order('left'); Ol
              Order of Quaternion Algebra (-1, -11) with base ring Rational Field with basis (1/2 + 1/2*j + 35*k, 1/4*i + 1/2*j + 75/4*k, j + 32*k, 38*k)
              sage: Or = I._compute_order('right'); Or
              Order of Quaternion Algebra (-1, -11) with base ring Rational Field with basis (1/2 + 1/2*j + 16*k, 1/2*i + 11/2*k, j + 13*k, 19*k)
              sage: Ol.discriminant()
              209
              sage: Or.discriminant()
              209
              sage: I.left_order() == Ol
              True
              sage: I.right_order() == Or
              True
  
          ALGORITHM: Let `b_1, b_2, b_3, b_3` be a basis for this
          fractional ideal `I`, and assume we want to compute the left
          order of `I` in the quaternion algebra `Q`.  Then
          multiplication by `b_i` on the right defines a map `B_i:Q \to
          Q`.  We have
  
          .. MATH::
  
             R = B_1^{-1}(I) \cap B_2^{-1}(I) \cap B_3^{-1}(I)\cap B_4^{-1}(I).
  
          This is because
  
          .. MATH::
  
             B_n^{-1}(I) = \{\alpha \in Q : \alpha b_n \in I \},
  
          and
  
          .. MATH::
  
             R = \{\alpha \in Q : \alpha b_n \in I, n=1,2,3,4\}.
          """
          if side == 'left':
              action = 'right'
          elif side == 'right':
              action = 'left'
          else:
              raise ValueError("side must be 'left' or 'right'")
          Q = self.quaternion_algebra()
          if Q.base_ring() != QQ:
              raise NotImplementedError("computation of left and right orders only implemented over QQ")
          M = [(~b).matrix(action=action) for b in self.basis()]
          B = self.basis_matrix()
          invs = [B * m for m in M]
          # Now intersect the row spans of each matrix in invs
          ISB = [Q(v) for v in intersection_of_row_modules_over_ZZ(invs).row_module(ZZ).basis()]
          return Q.quaternion_order(ISB)
  
      def left_order(self):
          """
          Return the left order associated to this fractional ideal.
  
          OUTPUT: an order in a quaternion algebra
  
          EXAMPLES::
  
              sage: B = BrandtModule(11)
              sage: R = B.maximal_order()
              sage: I = R.unit_ideal()
              sage: I.left_order()
              Order of Quaternion Algebra (-1, -11) with base ring Rational Field with basis (1/2 + 1/2*j, 1/2*i + 1/2*k, j, k)
  
          We do a consistency check::
  
              sage: B = BrandtModule(11,19); R = B.right_ideals()
              sage: [r.left_order().discriminant() for r in R]
              [209, 209, 209, 209, 209, 209, 209, 209, 209, 209, 209, 209, 209, 209, 209, 209, 209, 209]
          """
          if self.__left_order is None:
              self.__left_order = self._compute_order(side='left')
          return self.__left_order
  
      def right_order(self):
          """
          Return the right order associated to this fractional ideal.
  
          OUTPUT: an order in a quaternion algebra
  
          EXAMPLES::
  
              sage: I = BrandtModule(389).right_ideals()[1]; I
              Fractional ideal (2 + 6*j + 2*k, i + 2*j + k, 8*j, 8*k)
              sage: I.right_order()
              Order of Quaternion Algebra (-2, -389) with base ring Rational Field with basis (1/2 + 1/2*j + 1/2*k, 1/4*i + 1/2*j + 1/4*k, j, k)
              sage: I.left_order()
              Order of Quaternion Algebra (-2, -389) with base ring Rational Field with basis (1/2 + 1/2*j + 3/2*k, 1/8*i + 1/4*j + 9/8*k, j + k, 2*k)
  
          The following is a big consistency check.  We take reps for
          all the right ideal classes of a certain order, take the
          corresponding left orders, then take ideals in the left orders
          and from those compute the right order again::
  
              sage: B = BrandtModule(11,19); R = B.right_ideals()
              sage: O = [r.left_order() for r in R]
              sage: J = [O[i].left_ideal(R[i].basis()) for i in range(len(R))]
              sage: len(set(J))
              18
              sage: len(set([I.right_order() for I in J]))
              1
              sage: J[0].right_order() == B.order_of_level_N()
              True
          """
          if self.__right_order is None:
              self.__right_order = self._compute_order(side='right')
          return self.__right_order
  
      def __repr__(self):
          """
          Return string representation of this quaternion fractional ideal.
  
          EXAMPLES::
  
              sage: I = BrandtModule(11).right_ideals()[1]
              sage: type(I)
              <class 'sage.algebras.quatalg.quaternion_algebra.QuaternionFractionalIdeal_rational'>
              sage: I.__repr__()
              'Fractional ideal (2 + 6*j + 4*k, 2*i + 4*j + 2*k, 8*j, 8*k)'
          """
          return 'Fractional ideal %s' % (self.gens(),)
  
      def quaternion_order(self):
          """
          Return the order for which this ideal is a left or right
          fractional ideal.
  
          If this ideal has both a left and right
          ideal structure, then the left order is returned.  If it has
          neither structure, then an error is raised.
  
          OUTPUT: QuaternionOrder
  
          EXAMPLES::
  
              sage: R = QuaternionAlgebra(-11,-1).maximal_order()
              sage: R.unit_ideal().quaternion_order() is R
              doctest:...:  DeprecationWarning: quaternion_order() is deprecated, please use left_order() or right_order()
              See https://trac.sagemath.org/31583 for details.
              True
          """
          from sage.misc.superseded import deprecation
          deprecation(31583, 'quaternion_order() is deprecated, please use left_order() or right_order()')
          try:
              return self.__quaternion_order
          except AttributeError:
              pass
          if self.__left_order is not None:
              A = self.__left_order
          elif self.__right_order is not None:
              A = self.__right_order
          else:
              raise RuntimeError("unable to determine quaternion order of ideal without known order")
          self.__quaternion_order = A
          return A
  
      def ring(self):
          """
          Return ring that this is a fractional ideal for.
  
          The :meth:`ring` method will be removed from this class in the
          future.  Calling :meth:`ring` will then return the ambient
          quaternion algebra.  This is consistent with the behaviour for
          number fields.
  
          EXAMPLES::
  
              sage: R = QuaternionAlgebra(-11,-1).maximal_order()
              sage: R.unit_ideal().ring() is R
              doctest:...:  DeprecationWarning: ring() will return the quaternion algebra in the future, please use left_order() or right_order()
              See https://trac.sagemath.org/31583 for details.
              True
          """
          from sage.misc.superseded import deprecation
          deprecation(31583, 'ring() will return the quaternion algebra in the future, please use left_order() or right_order()')
          if self.__left_order is not None:
              return self.__left_order
          elif self.__right_order is not None:
              return self.__right_order
          else:
              raise RuntimeError("unable to determine quaternion order of ideal without known order")
  
      def basis(self):
          """
          Return a basis for this fractional ideal.
  
          OUTPUT: tuple
  
          EXAMPLES::
  
              sage: QuaternionAlgebra(-11,-1).maximal_order().unit_ideal().basis()
              (1/2 + 1/2*i, 1/2*j - 1/2*k, i, -k)
          """
          return self.gens()
  
      def _richcmp_(self, right, op):
          """
          Compare this fractional quaternion ideal to ``right``.
  
          EXAMPLES::
  
              sage: I = QuaternionAlgebra(-11, -1).maximal_order().unit_ideal()
              sage: I == I                # indirect doctest
              True
              sage: I == 5
              False
  
              sage: J = QuaternionAlgebra(-7, -1).maximal_order().unit_ideal()
              sage: J == I
              False
  
          Ideals can be equal even if they are defined by different
          bases (see :trac:`32245`)::
  
              sage: I == I.scale(-1)
              True
  
              sage: I != I                # indirect doctest
              False
          """
          return self.basis_matrix()._richcmp_(right.basis_matrix(), op)
  
      def __hash__(self):
          """
          Return the hash of ``self``.
  
          EXAMPLES::
  
              sage: I = QuaternionAlgebra(-11,-1).maximal_order().unit_ideal()
              sage: hash(I) == hash(I)
              True
  
          TESTS::
  
              sage: R = QuaternionAlgebra(-11,-1).maximal_order()
              sage: H = hash(R.right_ideal(R.basis()))
          """
          return hash(self.gens())
  
      @cached_method
      def basis_matrix(self):
          r"""
          Return basis matrix `M` in Hermite normal form for self as a
          matrix with rational entries.
  
          If `Q` is the ambient quaternion algebra, then the `\ZZ`-span of
          the rows of `M` viewed as linear combinations of Q.basis() =
          `[1,i,j,k]` is the fractional ideal self.  Also,
          ``M * M.denominator()`` is an integer matrix in Hermite normal form.
  
          OUTPUT: matrix over `\QQ`
  
          EXAMPLES::
  
              sage: QuaternionAlgebra(-11,-1).maximal_order().unit_ideal().basis_matrix()
              [1/2 1/2   0   0]
              [  0   1   0   0]
              [  0   0 1/2 1/2]
              [  0   0   0   1]
          """
          B = quaternion_algebra_cython.rational_matrix_from_rational_quaternions(self.gens())
          C, d = B._clear_denom()
          return C.hermite_form() / d
  
      def theta_series_vector(self, B):
          r"""
          Return theta series coefficients of ``self``, as a vector
          of ``B`` integers.
  
          INPUT:
  
          - ``B`` -- positive integer
  
          OUTPUT:
  
          Vector over `\ZZ` with ``B`` entries.
  
          EXAMPLES::
  
              sage: I = BrandtModule(37).right_ideals()[1]; I
              Fractional ideal (2 + 6*j + 2*k, i + 2*j + k, 8*j, 8*k)
              sage: I.theta_series_vector(5)
              (1, 0, 2, 2, 6)
              sage: I.theta_series_vector(10)
              (1, 0, 2, 2, 6, 4, 8, 6, 10, 10)
              sage: I.theta_series_vector(5)
              (1, 0, 2, 2, 6)
          """
          B = Integer(B)
          try:
              if len(self.__theta_series_vector) >= B:
                  return self.__theta_series_vector[:B]
          except AttributeError:
              pass
          V = FreeModule(ZZ, B)
          Q = self.quadratic_form()
          v = V(Q.representation_number_list(B))
          self.__theta_series_vector = v
          return v
  
      @cached_method
      def quadratic_form(self):
          """
          Return the normalized quadratic form associated to this quaternion ideal.
  
          OUTPUT: quadratic form
  
          EXAMPLES::
  
              sage: I = BrandtModule(11).right_ideals()[1]
              sage: Q = I.quadratic_form(); Q
              Quadratic form in 4 variables over Rational Field with coefficients:
              [ 18 22 33 22 ]
              [ * 7 22 11 ]
              [ * * 22 0 ]
              [ * * * 22 ]
              sage: Q.theta_series(10)
              1 + 12*q^2 + 12*q^3 + 12*q^4 + 12*q^5 + 24*q^6 + 24*q^7 + 36*q^8 + 36*q^9 + O(q^10)
              sage: I.theta_series(10)
              1 + 12*q^2 + 12*q^3 + 12*q^4 + 12*q^5 + 24*q^6 + 24*q^7 + 36*q^8 + 36*q^9 + O(q^10)
          """
          from sage.quadratic_forms.quadratic_form import QuadraticForm
          # first get the gram matrix
          gram_matrix = self.gram_matrix()
          # rescale so that there are no denominators
          gram_matrix, _ = gram_matrix._clear_denom()
          # Make sure gcd of all entries is 1.
          g = gram_matrix.gcd()
          if g != 1:
              gram_matrix = gram_matrix / g
          # now get the quadratic form
          return QuadraticForm(gram_matrix)
  
      def theta_series(self, B, var='q'):
          r"""
          Return normalized theta series of ``self``, as a power series over
          `\ZZ` in the variable ``var``, which is 'q' by default.
  
          The normalized theta series is by definition
  
          .. MATH::
  
              \theta_I(q) = \sum_{x \in I} q^{\frac{N(x)}{N(I)}}.
  
          INPUT:
  
          - ``B`` -- positive integer
          - ``var`` -- string (default: 'q')
  
          OUTPUT: power series
  
          EXAMPLES::
  
              sage: I = BrandtModule(11).right_ideals()[1]; I
              Fractional ideal (2 + 6*j + 4*k, 2*i + 4*j + 2*k, 8*j, 8*k)
              sage: I.norm()
              32
              sage: I.theta_series(5)
              1 + 12*q^2 + 12*q^3 + 12*q^4 + O(q^5)
              sage: I.theta_series(5,'T')
              1 + 12*T^2 + 12*T^3 + 12*T^4 + O(T^5)
              sage: I.theta_series(3)
              1 + 12*q^2 + O(q^3)
          """
          try:
              if self.__theta_series.prec() >= B:
                  if var == self.__theta_series.variable():
                      return self.__theta_series.add_bigoh(B)
                  else:
                      p_ring = self._theta_series.parent().change_var(var)
                      p_ring(self.__theta_series.list()[:B+1])
          except AttributeError:
              pass
          v = self.theta_series_vector(B)
          p_ring = PowerSeriesRing(ZZ, var)
          theta = p_ring(v.list()).add_bigoh(B)
          self.__theta_series = theta
          return theta
  
      @cached_method
      def gram_matrix(self):
          r"""
          Return the Gram matrix of this fractional ideal.
  
          OUTPUT: `4 \times 4` matrix over `\QQ`.
  
          EXAMPLES::
  
              sage: I = BrandtModule(3,5).right_ideals()[1]; I
              Fractional ideal (2 + 6*j + 4*k, 2*i + 4*j + 34*k, 8*j + 32*k, 40*k)
              sage: I.gram_matrix()
              [  640  1920  2112  1920]
              [ 1920 14080 13440 16320]
              [ 2112 13440 13056 15360]
              [ 1920 16320 15360 19200]
          """
          A = self.gens()
          two = QQ(2)
          m = [two * a.pair(b) for b in A for a in A]
          M44 = MatrixSpace(QQ, 4, 4)
          return M44(m, coerce=False)
  
      def norm(self):
          """
          Return the reduced norm of this fractional ideal.
  
          OUTPUT: rational number
  
          EXAMPLES::
  
              sage: M = BrandtModule(37)
              sage: C = M.right_ideals()
              sage: [I.norm() for I in C]
              [16, 32, 32]
  
              sage: (a,b) = M.quaternion_algebra().invariants()                                       # optional - magma
              sage: magma.eval('A<i,j,k> := QuaternionAlgebra<Rationals() | %s, %s>' % (a,b))         # optional - magma
              ''
              sage: magma.eval('O := QuaternionOrder(%s)' % str(list(C[0].right_order().basis())))    # optional - magma
              ''
              sage: [ magma('rideal<O | %s>' % str(list(I.basis()))).Norm() for I in C]               # optional - magma
              [16, 32, 32]
  
              sage: A.<i,j,k> = QuaternionAlgebra(-1,-1)
              sage: R = A.ideal([i,j,k,1/2 + 1/2*i + 1/2*j + 1/2*k])      # this is actually an order, so has reduced norm 1
              sage: R.norm()
              1
              sage: [ J.norm() for J in R.cyclic_right_subideals(3) ]     # enumerate maximal right R-ideals of reduced norm 3, verify their norms
              [3, 3, 3, 3]
          """
          G = self.gram_matrix() / QQ(2)
          r = G.det().abs()
          assert r.is_square(), "first is bad!"
          r = r.sqrt()
          # If we know either the left- or the right order, use that one to compute the norm.
          R = self.__left_order or self.__right_order or self.left_order()
          r /= R.discriminant()
          assert r.is_square(), "second is bad!"
          return r.sqrt()
  
      def conjugate(self):
          """
          Return the ideal with generators the conjugates of the generators for self.
  
          OUTPUT: a quaternionic fractional ideal
  
          EXAMPLES::
  
              sage: I = BrandtModule(3,5).right_ideals()[1]; I
              Fractional ideal (2 + 6*j + 4*k, 2*i + 4*j + 34*k, 8*j + 32*k, 40*k)
              sage: I.conjugate()
              Fractional ideal (2 + 2*j + 28*k, 2*i + 4*j + 34*k, 8*j + 32*k, 40*k)
          """
          return self.quaternion_algebra().ideal([b.conjugate() for b in self.basis()],
                                                 left_order=self.__right_order,
                                                 right_order=self.__left_order)
  
      def __mul__(self, right):
          """
          Return the product of the fractional ideals ``self`` and ``right``.
  
          .. note::
  
             We do not keep track of left or right order structure.
  
          EXAMPLES::
  
              sage: I = BrandtModule(3,5).right_ideals()[1]; I
              Fractional ideal (2 + 6*j + 4*k, 2*i + 4*j + 34*k, 8*j + 32*k, 40*k)
              sage: I*I
              Fractional ideal (8 + 24*j + 16*k, 8*i + 16*j + 136*k, 32*j + 128*k, 160*k)
              sage: I*I.conjugate()
              Fractional ideal (16 + 16*j + 224*k, 8*i + 16*j + 136*k, 32*j + 128*k, 320*k)
              sage: I.multiply_by_conjugate(I)
              Fractional ideal (16 + 16*j + 224*k, 8*i + 16*j + 136*k, 32*j + 128*k, 320*k)
          """
          if isinstance(right, QuaternionOrder):
              right = right.unit_ideal()
          if not isinstance(right, QuaternionFractionalIdeal_rational):
              return self.scale(right, left=False)
          gens = [a*b for a in self.basis() for b in right.basis()]
          # if self.__right_order == right.__left_order:
          #     left_order = self.__left_order
          #     right_order = right.__right_order
          basis = tuple(basis_for_quaternion_lattice(gens))
          A = self.quaternion_algebra()
          return A.ideal(basis, check=False)
  
      def __add__(self, other):
          """
          Return the sum of the fractional ideals ``self`` and ``other``.
  
          EXAMPLES::
  
              sage: I = BrandtModule(11,5).right_ideals()[1]; I
              Fractional ideal (2 + 2*j + 20*k, 2*i + 4*j + 6*k, 8*j, 40*k)
              sage: J = BrandtModule(11,5).right_ideals()[2]; J
              Fractional ideal (2 + 6*j + 20*k, 2*i + 4*j + 26*k, 8*j, 40*k)
              sage: I + J
              Fractional ideal (2 + 2*j, 2*i + 6*k, 4*j, 20*k)
          """
          if isinstance(other, QuaternionOrder):
              other = other.unit_ideal()
          if not isinstance(other, QuaternionFractionalIdeal_rational):
              raise TypeError("can only add quaternion ideals")
          return self.quaternion_algebra().ideal(self.basis() + other.basis())
  
      def _acted_upon_(self, other, on_left):
          """
          Scale a quaternion ideal.
  
          EXAMPLES::
  
              sage: Q.<i,j,k> = QuaternionAlgebra(-1,-419)
              sage: O = Q.maximal_order()
              sage: I = 7*O.unit_ideal() + O.unit_ideal()*(j-1); I
              Fractional ideal (1/2 + 13/2*j, 1/2*i + 13/2*k, 7*j, 7*k)
  
          TESTS::
  
              sage: (5+i-j)*I == I.scale(5+i-j, left=True)
              True
              sage: I*(5+i-j) == I.scale(5+i-j, left=False)
              True
          """
          return self.scale(other, left=(not on_left))
  
      @cached_method
      def free_module(self):
          r"""
          Return the underlying free `\ZZ`-module corresponding to this ideal.
  
          OUTPUT:
  
          Free `\ZZ`-module of rank 4 embedded in an ambient `\QQ^4`.
  
          EXAMPLES::
  
              sage: X = BrandtModule(3,5).right_ideals()
              sage: X[0]
              Fractional ideal (2 + 2*j + 8*k, 2*i + 18*k, 4*j + 16*k, 20*k)
              sage: X[0].free_module()
              Free module of degree 4 and rank 4 over Integer Ring
              Echelon basis matrix:
              [ 2  0  2  8]
              [ 0  2  0 18]
              [ 0  0  4 16]
              [ 0  0  0 20]
              sage: X[0].scale(1/7).free_module()
              Free module of degree 4 and rank 4 over Integer Ring
              Echelon basis matrix:
              [ 2/7    0  2/7  8/7]
              [   0  2/7    0 18/7]
              [   0    0  4/7 16/7]
              [   0    0    0 20/7]
  
              sage: QuaternionAlgebra(-11,-1).maximal_order().unit_ideal().basis_matrix()
              [1/2 1/2   0   0]
              [  0   1   0   0]
              [  0   0 1/2 1/2]
              [  0   0   0   1]
  
          The free module method is also useful since it allows for checking if
          one ideal is contained in another, computing quotients `I/J`, etc.::
  
              sage: X = BrandtModule(3,17).right_ideals()
              sage: I = X[0].intersection(X[2]); I
              Fractional ideal (2 + 2*j + 164*k, 2*i + 4*j + 46*k, 16*j + 224*k, 272*k)
              sage: I.free_module().is_submodule(X[3].free_module())
              False
              sage: I.free_module().is_submodule(X[1].free_module())
              True
              sage: X[0].free_module() / I.free_module()
              Finitely generated module V/W over Integer Ring with invariants (4, 4)
  
          This shows that the issue at :trac:`6760` is fixed::
  
              sage: R.<i,j,k> = QuaternionAlgebra(-1, -13)
              sage: I = R.ideal([2+i, 3*i, 5*j, j+k]); I
              Fractional ideal (2 + i, 3*i, j + k, 5*k)
              sage: I.free_module()
              Free module of degree 4 and rank 4 over Integer Ring
              Echelon basis matrix:
              [2 1 0 0]
              [0 3 0 0]
              [0 0 1 1]
              [0 0 0 5]
          """
          return self.basis_matrix().row_module(ZZ)
  
      def intersection(self, J):
          """
          Return the intersection of the ideals self and `J`.
  
          EXAMPLES::
  
              sage: X = BrandtModule(3,5).right_ideals()
              sage: I = X[0].intersection(X[1]); I
              Fractional ideal (2 + 6*j + 4*k, 2*i + 4*j + 34*k, 8*j + 32*k, 40*k)
  
          """
          V = self.free_module().intersection(J.free_module())
          H, d = V.basis_matrix()._clear_denom()
          A = self.quaternion_algebra()
          gens = quaternion_algebra_cython.rational_quaternions_from_integral_matrix_and_denom(A, H, d)
          return A.ideal(gens)
  
      def multiply_by_conjugate(self, J):
          """
          Return product of self and the conjugate Jbar of `J`.
  
          INPUT:
  
          - ``J`` -- a quaternion ideal.
  
          OUTPUT: a quaternionic fractional ideal.
  
          EXAMPLES::
  
              sage: R = BrandtModule(3,5).right_ideals()
              sage: R[0].multiply_by_conjugate(R[1])
              Fractional ideal (8 + 8*j + 112*k, 8*i + 16*j + 136*k, 32*j + 128*k, 160*k)
              sage: R[0]*R[1].conjugate()
              Fractional ideal (8 + 8*j + 112*k, 8*i + 16*j + 136*k, 32*j + 128*k, 160*k)
          """
          Jbar = [b.conjugate() for b in J.basis()]
          gens = [a * b for a in self.basis() for b in Jbar]
          basis = tuple(basis_for_quaternion_lattice(gens))
          R = self.quaternion_algebra()
          return R.ideal(basis, check=False)
  
      def is_equivalent(I, J, B=10):
          """
          Return ``True`` if ``I`` and ``J`` are equivalent as right ideals.
  
          INPUT:
  
          - ``I`` -- a fractional quaternion ideal (self)
  
          - ``J`` -- a fractional quaternion ideal with same order as ``I``
  
          - ``B`` -- a bound to compute and compare theta series before
            doing the full equivalence test
  
          OUTPUT: bool
  
          EXAMPLES::
  
              sage: R = BrandtModule(3,5).right_ideals(); len(R)
              2
              sage: R[0].is_equivalent(R[1])
              False
              sage: R[0].is_equivalent(R[0])
              True
              sage: OO = R[0].left_order()
              sage: S = OO.right_ideal([3*a for a in R[0].basis()])
              sage: R[0].is_equivalent(S)
              True
          """
          if not isinstance(I, QuaternionFractionalIdeal_rational):
              return False
  
          if I.right_order() != J.right_order():
              raise ValueError("I and J must be right ideals")
  
          # Just test theta series first.  If the theta series are
          # different, the ideals are definitely not equivalent.
          if B > 0 and I.theta_series_vector(B) != J.theta_series_vector(B):
              return False
  
          # The theta series are the same, so perhaps the ideals are
          # equivalent.  We use Prop 1.18 of [Pizer, 1980] to decide.
          # 1. Compute I * Jbar
          # see Prop. 1.17 in Pizer.  Note that we use IJbar instead of
          # JbarI since we work with right ideals
          IJbar = I.multiply_by_conjugate(J)
  
          # 2. Determine if there is alpha in K such
          #    that N(alpha) = N(I)*N(J) as explained by Pizer.
          c = IJbar.theta_series_vector(2)[1]
          return c != 0
  
      def __contains__(self, x):
          """
          Return whether x is in self.
  
          EXAMPLES::
  
              sage: R.<i,j,k> = QuaternionAlgebra(-3, -13)
              sage: I = R.ideal([2+i, 3*i, 5*j, j+k])
              sage: 2+i in I
              True
              sage: 2+i+j+k in I
              True
              sage: 1+i in I
              False
              sage: 101*j + k in I
              True
          """
          try:
              x = self.quaternion_algebra()(x)
              return self.basis_matrix().transpose().solve_right(vector(x)) in ZZ**4
          except (ValueError, TypeError):
              return False
  
      @cached_method
      def cyclic_right_subideals(self, p, alpha=None):
          r"""
          Let `I` = ``self``. This function returns the right subideals
          `J` of `I` such that `I/J` is an `\GF{p}`-vector space of
          dimension 2.
  
          INPUT:
  
          - ``p`` -- prime number (see below)
  
          - ``alpha`` -- (default: ``None``) element of quaternion algebra,
            which can be used to parameterize the order of the
            ideals `J`.  More precisely the `J`'s are the right annihilators
            of `(1,0) \alpha^i` for `i=0,1,2,...,p`
  
          OUTPUT:
  
          - list of right ideals
  
          .. NOTE::
  
             Currently, `p` must satisfy a bunch of conditions, or a
             ``NotImplementedError`` is raised.  In particular, `p` must be
             odd and unramified in the quaternion algebra, must be
             coprime to the index of the right order in the maximal
             order, and also coprime to the normal of self.  (The Brandt
             modules code has a more general algorithm in some cases.)
  
          EXAMPLES::
  
              sage: B = BrandtModule(2,37); I = B.right_ideals()[0]
              sage: I.cyclic_right_subideals(3)
              [Fractional ideal (2 + 2*i + 10*j + 90*k, 4*i + 4*j + 152*k, 12*j + 132*k, 444*k), Fractional ideal (2 + 2*i + 2*j + 150*k, 4*i + 8*j + 196*k, 12*j + 132*k, 444*k), Fractional ideal (2 + 2*i + 6*j + 194*k, 4*i + 8*j + 344*k, 12*j + 132*k, 444*k), Fractional ideal (2 + 2*i + 6*j + 46*k, 4*i + 4*j + 4*k, 12*j + 132*k, 444*k)]
  
              sage: B = BrandtModule(5,389); I = B.right_ideals()[0]
              sage: C = I.cyclic_right_subideals(3); C
              [Fractional ideal (2 + 10*j + 546*k, i + 6*j + 133*k, 12*j + 3456*k, 4668*k), Fractional ideal (2 + 2*j + 2910*k, i + 6*j + 3245*k, 12*j + 3456*k, 4668*k), Fractional ideal (2 + i + 2295*k, 3*i + 2*j + 3571*k, 4*j + 2708*k, 4668*k), Fractional ideal (2 + 2*i + 2*j + 4388*k, 3*i + 2*j + 2015*k, 4*j + 4264*k, 4668*k)]
              sage: [(I.free_module()/J.free_module()).invariants() for J in C]
              [(3, 3), (3, 3), (3, 3), (3, 3)]
              sage: I.scale(3).cyclic_right_subideals(3)
              [Fractional ideal (6 + 30*j + 1638*k, 3*i + 18*j + 399*k, 36*j + 10368*k, 14004*k), Fractional ideal (6 + 6*j + 8730*k, 3*i + 18*j + 9735*k, 36*j + 10368*k, 14004*k), Fractional ideal (6 + 3*i + 6885*k, 9*i + 6*j + 10713*k, 12*j + 8124*k, 14004*k), Fractional ideal (6 + 6*i + 6*j + 13164*k, 9*i + 6*j + 6045*k, 12*j + 12792*k, 14004*k)]
              sage: C = I.scale(1/9).cyclic_right_subideals(3); C
              [Fractional ideal (2/9 + 10/9*j + 182/3*k, 1/9*i + 2/3*j + 133/9*k, 4/3*j + 384*k, 1556/3*k), Fractional ideal (2/9 + 2/9*j + 970/3*k, 1/9*i + 2/3*j + 3245/9*k, 4/3*j + 384*k, 1556/3*k), Fractional ideal (2/9 + 1/9*i + 255*k, 1/3*i + 2/9*j + 3571/9*k, 4/9*j + 2708/9*k, 1556/3*k), Fractional ideal (2/9 + 2/9*i + 2/9*j + 4388/9*k, 1/3*i + 2/9*j + 2015/9*k, 4/9*j + 4264/9*k, 1556/3*k)]
              sage: [(I.scale(1/9).free_module()/J.free_module()).invariants() for J in C]
              [(3, 3), (3, 3), (3, 3), (3, 3)]
  
              sage: Q.<i,j,k> = QuaternionAlgebra(-2,-5)
              sage: I = Q.ideal([Q(1),i,j,k])
              sage: I.cyclic_right_subideals(3)
              [Fractional ideal (1 + 2*j, i + k, 3*j, 3*k), Fractional ideal (1 + j, i + 2*k, 3*j, 3*k), Fractional ideal (1 + 2*i, 3*i, j + 2*k, 3*k), Fractional ideal (1 + i, 3*i, j + k, 3*k)]
  
          The general algorithm is not yet implemented here::
  
              sage: I.cyclic_right_subideals(3)[0].cyclic_right_subideals(3)
              Traceback (most recent call last):
              ...
              NotImplementedError: general algorithm not implemented (The given basis vectors must be linearly independent.)
          """
          R = self.right_order()
          Q = self.quaternion_algebra()
          f = Q.modp_splitting_map(p)
          if alpha is not None:
              alpha = f(alpha)
          W = GF(p)**4
          try:
              A = W.span_of_basis([W(f(a).list()) for a in self.basis()])
              scale = 1
              IB = self.basis_matrix()
          except (ValueError, ZeroDivisionError):
              # try rescaling the ideal.
              B, d = self.basis_matrix()._clear_denom()
              g = gcd(B.list())
              IB = B / g
              scale = g / d
              try:
                  A = W.span_of_basis([W(f(Q(a.list())).list()) for a in IB.rows()])
              except (ValueError, ZeroDivisionError) as msg:
                  # Here we could replace the ideal by an *equivalent*
                  # ideal that works.  This is always possible.
                  # However, I haven't implemented that algorithm yet.
                  raise NotImplementedError("general algorithm not implemented (%s)" % msg)
  
          Ai = A.basis_matrix()**(-1)
          AiB = Ai.change_ring(QQ) * IB
  
          # Do not care about the denominator since we're really working in I/p*I.
          AiB, _ = AiB._clear_denom()
  
          pB = p*IB
          pB, d = pB._clear_denom()
  
          ans = []
          Z = matrix(ZZ, 2, 4)
  
          P1 = P1List(p)
          if alpha is None:
              lines = P1
          else:
              x = alpha
              lines = []
              for _ in range(p + 1):
                  lines.append(P1.normalize(x[0, 0], x[0, 1]))
                  x *= alpha
  
          for u, v in lines:
              # The following does:
              #    z = matrix(QQ,2,4,[0,-v,0,u, -v,0,u,0],check=False) * AiB
              Z[0, 1] = -v
              Z[0, 3] = u
              Z[1, 0] = -v
              Z[1, 2] = u
              z = Z * AiB
              # Now construct submodule of the ideal I spanned by the
              # linear combinations given by z of the basis for J along
              # with p*I.
              G = (d*z).stack(pB)  # have to multiply by d since we divide by it below in the "gens = " line.
              H = G._hnf_pari(0, include_zero_rows=False)
              gens = tuple(quaternion_algebra_cython.rational_quaternions_from_integral_matrix_and_denom(Q, H, d))
              if scale != 1:
                  gens = tuple([scale * gg for gg in gens])
              J = R.right_ideal(gens, check=False)
              ans.append(J)
          return ans
  
  #######################################################################
  # Some utility functions that are needed here and are too
  # specialized to go elsewhere.
  #######################################################################
  
  
  def basis_for_quaternion_lattice(gens, reverse=False):
      r"""
      Return a basis for the `\ZZ`-lattice in a quaternion algebra
      spanned by the given gens.
  
      INPUT:
  
      - ``gens`` -- list of elements of a single quaternion algebra
  
      - ``reverse`` -- when computing the HNF do it on the basis
        `(k,j,i,1)` instead of `(1,i,j,k)`; this ensures
        that if ``gens`` are the generators for an order,
        the first returned basis vector is 1
  
      EXAMPLES::
  
          sage: from sage.algebras.quatalg.quaternion_algebra import basis_for_quaternion_lattice
          sage: A.<i,j,k> = QuaternionAlgebra(-1,-7)
          sage: basis_for_quaternion_lattice([i+j, i-j, 2*k, A(1/3)])
          [1/3, i + j, 2*j, 2*k]
  
          sage: basis_for_quaternion_lattice([A(1),i,j,k])
          [1, i, j, k]
  
      """
      if not gens:
          return []
      Z, d = quaternion_algebra_cython.integral_matrix_and_denom_from_rational_quaternions(gens, reverse)
      H = Z._hnf_pari(0, include_zero_rows=False)
      A = gens[0].parent()
      return quaternion_algebra_cython.rational_quaternions_from_integral_matrix_and_denom(A, H, d, reverse)
  
  
  def intersection_of_row_modules_over_ZZ(v):
      r"""
      Intersect the `\ZZ`-modules with basis matrices the full rank `4 \times 4`
      `\QQ`-matrices in the list v.
  
      The returned intersection is
      represented by a `4 \times 4` matrix over `\QQ`.  This can also be done
      using modules and intersection, but that would take over twice as long
      because of overhead, hence this function.
  
      EXAMPLES::
  
          sage: a = matrix(QQ,4,[-2, 0, 0, 0, 0, -1, -1, 1, 2, -1/2, 0, 0, 1, 1, -1, 0])
          sage: b = matrix(QQ,4,[0, -1/2, 0, -1/2, 2, 1/2, -1, -1/2, 1, 2, 1, -2, 0, -1/2, -2, 0])
          sage: c = matrix(QQ,4,[0, 1, 0, -1/2, 0, 0, 2, 2, 0, -1/2, 1/2, -1, 1, -1, -1/2, 0])
          sage: v = [a,b,c]
          sage: from sage.algebras.quatalg.quaternion_algebra import intersection_of_row_modules_over_ZZ
          sage: M = intersection_of_row_modules_over_ZZ(v); M
          [   2    0   -1   -1]
          [  -4    1    1   -3]
          [  -3 19/2   -1   -4]
          [   2   -3   -8    4]
          sage: M2 = a.row_module(ZZ).intersection(b.row_module(ZZ)).intersection(c.row_module(ZZ))
          sage: M.row_module(ZZ) == M2
          True
      """
      if len(v) <= 0:
          raise ValueError("v must have positive length")
      if len(v) == 1:
          return v[0]
      elif len(v) == 2:
          # real work - the base case
          a, b = v
          s, _ = a.stack(b)._clear_denom()
          s = s.transpose()
          K = s.right_kernel_matrix(algorithm='pari', basis='computed')
          n = a.nrows()
          return K.matrix_from_columns(range(n)) * a
      else:
          # induct
          w = intersection_of_row_modules_over_ZZ(v[:2])
          return intersection_of_row_modules_over_ZZ([w] + v[2:])
  
  
  def normalize_basis_at_p(e, p, B=QuaternionAlgebraElement_abstract.pair):
      r"""
      Compute a (at ``p``) normalized basis from the given basis ``e``
      of a `\ZZ`-module.
  
      The returned basis is (at ``p``) a `\ZZ_p` basis for the same
      module, and has the property that with respect to it the quadratic
      form induced by the bilinear form B is represented as a orthogonal
      sum of atomic forms multiplied by p-powers.
  
      If `p \neq 2` this means that the form is diagonal with respect to
      this basis.
  
      If `p = 2` there may be additional 2-dimensional subspaces on which
      the form is represented as `2^e (ax^2 + bxy + cx^2)` with
      `0 = v_2(b) = v_2(a) \leq v_2(c)`.
  
      INPUT:
  
      - ``e`` -- list; basis of a `\ZZ` module.
        WARNING: will be modified!
  
      - ``p`` -- prime for at which the basis should be normalized
  
      - ``B`` --
        (default: :meth:`QuaternionAlgebraElement_abstract.pair`)
        a bilinear form with respect to which to normalize
  
      OUTPUT:
  
      - A list containing two-element tuples: The first element of
        each tuple is a basis element, the second the valuation of
        the orthogonal summand to which it belongs. The list is sorted
        by ascending valuation.
  
      EXAMPLES::
  
          sage: from sage.algebras.quatalg.quaternion_algebra import normalize_basis_at_p
          sage: A.<i,j,k> = QuaternionAlgebra(-1, -1)
          sage: e = [A(1), i, j, k]
          sage: normalize_basis_at_p(e, 2)
          [(1, 0), (i, 0), (j, 0), (k, 0)]
  
          sage: A.<i,j,k> = QuaternionAlgebra(210)
          sage: e = [A(1), i, j, k]
          sage: normalize_basis_at_p(e, 2)
          [(1, 0), (i, 1), (j, 1), (k, 2)]
  
          sage: A.<i,j,k> = QuaternionAlgebra(286)
          sage: e = [A(1), k, 1/2*j + 1/2*k, 1/2 + 1/2*i + 1/2*k]
          sage: normalize_basis_at_p(e, 5)
          [(1, 0), (1/2*j + 1/2*k, 0), (-5/6*j + 1/6*k, 1), (1/2*i, 1)]
  
          sage: A.<i,j,k> = QuaternionAlgebra(-1,-7)
          sage: e = [A(1), k, j, 1/2 + 1/2*i + 1/2*j + 1/2*k]
          sage: normalize_basis_at_p(e, 2)
          [(1, 0), (1/2 + 1/2*i + 1/2*j + 1/2*k, 0), (-34/105*i - 463/735*j + 71/105*k, 1), (-34/105*i - 463/735*j + 71/105*k, 1)]
      """
  
      N = len(e)
      if N == 0:
          return []
      else:
          min_m, min_n, min_v = 0, 0, infinity
  
          # Find two basis vector on which the bilinear form has minimal
          # p-valuation. If there is more than one such pair, always
          # prefer diagonal entries over any other and (secondary) take
          # min_m and then min_n as small as possible
          for m in range(N):
              for n in range(m, N):
                  v = B(e[m], e[n]).valuation(p)
                  if v < min_v or (v == min_v and (min_m != min_n) and (m == n)):
                      min_m, min_n, min_v = m, n, v
  
          if (min_m == min_n) or p != 2:      # In this case we can diagonalize
              if min_m == min_n:              # Diagonal entry has minimal valuation
                  f0 = e[min_m]
              else:
                  f0 = e[min_m] + e[min_n]    # Only off-diagonal entries have min. val., but p!=2
  
              # Swap with first vector
              e[0], e[min_m] = e[min_m], e[0]
  
              # Orthogonalize remaining vectors with respect to f
              c = B(f0, f0)
              for l in range(1, N):
                  e[l] = e[l] - B(e[l], f0) / c * f0
  
              # Recursively normalize remaining vectors
              f = normalize_basis_at_p(e[1:], p)
              f.insert(0, (f0, min_v - valuation(p, 2)))
              return f
  
          else:   # p = 2 and only off-diagonal entries have min. val., gives 2-dim. block
              # first diagonal entry should have smaller valuation
              if B(e[min_m], e[min_m]).valuation(p) > B(e[min_n], e[min_n]).valuation(p):
                  e[min_m], e[min_n] = e[min_n], e[min_m]
  
              f0 = p**min_v / B(e[min_m], e[min_n]) * e[min_m]
              f1 = e[min_n]
  
              # Ensures that (B(f0,f0)/2).valuation(p) <= B(f0,f1).valuation(p)
              if B(f0, f1).valuation(p) + 1 < B(f0, f0).valuation(p):
                  f0 += f1
                  f1 = f0
  
              # Make remaining vectors orthogonal to span of f0, f1
              e[min_m] = e[0]
              e[min_n] = e[1]
  
              B00 = B(f0, f0)
              B11 = B(f1, f1)
              B01 = B(f0, f1)
              d = B00 * B11 - B01**2
              tu = [(B01 * B(f1, e[l]) - B11 * B(f0, e[l]),
                     B01 * B(f0, e[l]) - B00 * B(f1, e[l])) for l in range(2, N)]
  
              e[2:n] = [e[l] + tu[l-2][0]/d * f0 + tu[l-2][1]/d * f1 for l in range(2, N)]
  
              # Recursively normalize remaining vectors
              f = normalize_basis_at_p(e[2:N], p)
              return [(f0, min_v), (f1, min_v)] + f
  
  
  def maxord_solve_aux_eq(a, b, p):
      r"""
      Given ``a`` and ``b`` and an even prime ideal ``p`` find
      (y,z,w) with y a unit mod `p^{2e}` such that
  
      .. MATH::
  
          1 - ay^2 - bz^2 + abw^2 \equiv 0 mod p^{2e},
  
      where `e` is the ramification index of `p`.
  
      Currently only `p=2` is implemented by hardcoding solutions.
  
      INPUT:
  
      - ``a`` -- integer with `v_p(a) = 0`
  
      - ``b`` -- integer with `v_p(b) \in \{0,1\}`
  
      - ``p`` -- even prime ideal (actually only ``p=ZZ(2)`` is implemented)
  
      OUTPUT:
  
      - A tuple `(y, z, w)`
  
      EXAMPLES::
  
          sage: from sage.algebras.quatalg.quaternion_algebra import maxord_solve_aux_eq
          sage: for a in [1,3]:
          ....:     for b in [1,2,3]:
          ....:         (y,z,w) = maxord_solve_aux_eq(a, b, 2)
          ....:         assert mod(y, 4) == 1 or mod(y, 4) == 3
          ....:         assert mod(1 - a*y^2 - b*z^2 + a*b*w^2, 4) == 0
      """
      if p != ZZ(2):
          raise NotImplementedError("Algorithm only implemented over ZZ at the moment")
  
      v_a = a.valuation(p)
      v_b = b.valuation(p)
  
      if v_a != 0:
          raise RuntimeError("a must have v_p(a)=0")
      if v_b != 0 and v_b != 1:
          raise RuntimeError("b must have v_p(b) in {0,1}")
  
      R = ZZ.quo(ZZ(4))
      lut = {(R(1), R(1)): (1, 1, 1),
             (R(1), R(2)): (1, 0, 0),
             (R(1), R(3)): (1, 0, 0),
             (R(3), R(1)): (1, 1, 1),
             (R(3), R(2)): (1, 0, 1),
             (R(3), R(3)): (1, 1, 1)}
  
      return lut[(R(a), R(b))]
  