Compute PARI permutation and other fixes

- Slightly adapted docstring for construction via ramification
- Added (currently non-deterministic and indirect) tests for correctness of PARI transfer, including "direction" of permutation
- Corrected order of constructor type checks
- Implemented permutation computation to give correct quaternion algebra with correct ramification (credit to @videlec for the elegant permutation method)

Amend:
- Avoid test failures due to non-deterministic calculations

src/sage/algebras/quatalg/quaternion_algebra.py

@@ -56,7 +56,9 @@
 from sage.rings.rational_field import is_RationalField, QQ
 from sage.rings.infinity import infinity
 from sage.rings.number_field.number_field_base import NumberField
+from sage.rings.qqbar import AA
 from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
+from sage.rings.polynomial.polynomial_ring import polygen
 from sage.rings.power_series_ring import PowerSeriesRing
 from sage.structure.category_object import normalize_names
 from sage.structure.parent import Parent
@@ -80,12 +82,12 @@
 from sage.modular.modsym.p1list import P1List
 
 from sage.misc.cachefunc import cached_method
-from sage.misc.functional import is_odd, category
+from sage.misc.functional import is_odd
 
 from sage.categories.algebras import Algebras
 from sage.categories.number_fields import NumberFields
 
-from sage.libs.pari.all import pari
+from sage.combinat.words.word import Word
 
 ########################################################
 # Constructor
@@ -111,12 +113,13 @@ class QuaternionAlgebraFactory(UniqueFactory):
       `D` over `K = \QQ`.  Suitable nonzero rational numbers `a`, `b`
       as above are deduced from `D`.
 
-    - ``QuaternionAlgebra(K, primes, inv_archimedean)``, where `primes`
-      is a list of prime ideals and `inv_archimedean` is a list of local
-      invariants (0 or 1/2) specifying the ramification at the (infinite)
-      real places of `K`.  This constructs a quaternion algebra ramified
-      exacly at the places in `primes` and those in `K.real_embeddings()`
-      indexed by `i` with `inv_archimedean[i]=1/2`.
+    - ``QuaternionAlgebra(K, primes, inv_archimedean)``, where `K` is a
+      number field or `\QQ`, `primes` is a list of prime ideals of `K`
+      and `inv_archimedean` is a list of local invariants (0 or 1/2)
+      specifying the ramification at the (infinite) real places of `K`.
+      This constructs a quaternion algebra ramified exacly at the places
+      given by `primes` and those in `K.embeddings(AA)` indexed by `i`
+      with `inv_archimedean[i]=1/2`.
 
     OUTPUT:
 
@@ -203,6 +206,8 @@ class QuaternionAlgebraFactory(UniqueFactory):
         Traceback (most recent call last):
         ...
         ValueError: quaternion algebra over the rationals must have an even number of ramified places
+
+        sage: x = polygen(ZZ, 'x')
         sage: K.<w> = NumberField(x^2-x-1)
         sage: P = K.prime_above(2)
         sage: Q = K.prime_above(3)
@@ -255,6 +260,28 @@ class QuaternionAlgebraFactory(UniqueFactory):
         Rational Field
         sage: parent(Q._b)
         Rational Field
+
+    Check that construction via ramification yields the correct algebra, i.e.
+    that the differences between Sage and PARI are incorporated correctly::
+        sage: x = polygen(ZZ, 'x')
+        sage: K.<v> = NumberField(-3*x^5 - 11*x^4 - 4*x^3 + 1)
+        sage: QuaternionAlgebra(K, [], [1/2, 0, 1/2])       # not tested
+        Quaternion Algebra (-87*v^4 - 412*v^3 - 472*v^2 - 173*v - 7, 6*v^4 + 16*v^3 - 26*v^2 - 40*v - 1)
+        with base ring Number Field in v with defining polynomial -3*x^5 - 11*x^4 - 4*x^3 + 1
+
+    Also check that the Sage-PARI permutation is the correct way around, i.e.
+    it does not need to be replaced by its inverse (computation is not deterministic,
+    this will be checked properly once the required code has been merged)::
+
+        sage: x = polygen(ZZ, 'x')
+        sage: K.<j> = NumberField(5*x^4 - 50*x^2 + 5)
+        sage: P = K.prime_above(2)
+        sage: Q = K.prime_above(5)
+        sage: inv_arch = [1/2, 1/2, 0, 0]
+        sage: QuaternionAlgebra(K, [P,Q], inv_arch)         # not tested
+        Quaternion Algebra (-51/4*j^3 - 11/4*j^2 + 523/4*j + 47/4, 149/2*j^3 + 85*j^2 - 1073/2*j - 900)
+        with base ring Number Field in j with defining polynomial 5*x^4 - 50*x^2 + 5            
+
     """
     def create_key(self, arg0, arg1=None, arg2=None, names='i,j,k'):
         """
@@ -293,37 +320,65 @@ def create_key(self, arg0, arg1=None, arg2=None, names='i,j,k'):
             a = K(v[0])
             b = K(v[1])
 
-        else:
+        elif isinstance(arg1, list) and isinstance(arg2, list):
             # QuaternionAlgebra(K, primes, inv_archimedean)
             K = arg0
-            if category(K) is not NumberFields():
-                raise ValueError("quaternion algbera must be defined over a number field")
-            if isinstance(arg1, list) and isinstance(arg2, list):
-                if not set(arg2).issubset(set([0,QQ(1/2)])):
-                    raise ValueError("list of local invariants specifying ramification should contain only 0 and 1/2")
-                arg1 = set(arg1)
-                if not all([p.is_prime() for p in arg1]):
-                    raise ValueError("quaternion algebra constructor requires a list of primes specifying the ramification")
-                if is_RationalField(K):
-                    if len(arg2) > 1 or (len(arg2) == 1 and is_odd(len(arg1) + 2*arg2[0])):
-                        raise ValueError("quaternion algebra over the rationals must have an even number of ramified places")
-                    D = ZZ.ideal_monoid().prod(arg1).gen()
-                    a, b = hilbert_conductor_inverse(D)
-                    a = QQ(a)
-                    b = QQ(b)
-                else:
-                    if len(arg2) != len(K.real_places()):
-                        raise ValueError("must specify ramification at the real places of %s" % K)
-                    if is_odd(len(arg1) + 2 * sum(arg2)):
-                        raise ValueError("quaternion algebra over a number field must have an even number of ramified places")
-                    fin_places_pari = [I.pari_prime() for I in arg1]
-                    A = pari(K).alginit([2, [fin_places_pari, [QQ(1/2)] * len(fin_places_pari)], arg2], maxord=0)
-                    a = K(A.algsplittingfield().disc()[1])
-                    b = K(A.algb())
+            if K not in NumberFields():
+                raise ValueError("quaternion algebra must be defined over a number field")
+            if not set(arg2).issubset(set([0, QQ(1/2)])):
+                raise ValueError("list of local invariants specifying ramification should contain only 0 and 1/2")
+            primes = set(arg1)
+            if not all([p.is_prime() for p in primes]):
+                raise ValueError("quaternion algebra constructor requires a list of primes specifying the ramification")
+            if is_RationalField(K):
+                if len(arg2) > 1 or (len(arg2) == 1 and is_odd(len(primes) + 2*arg2[0])):
+                    raise ValueError("quaternion algebra over the rationals must have an even number of ramified places")
+                D = ZZ.ideal_monoid().prod(primes).gen()
+                a, b = hilbert_conductor_inverse(D)
+                a = QQ(a)
+                b = QQ(b)
             else:
-                # QuaternionAlgebra(K, a, b)
-                a = K(arg1)
-                b = K(arg2)
+                if len(arg2) != len(K.real_places()):
+                    raise ValueError("must specify ramification at the real places of %s" % K)
+                if is_odd(len(primes) + 2 * sum(arg2)):
+                    raise ValueError("quaternion algebra over a number field must have an even number of ramified places")
+
+                # We want to compute the correct quaternion algebra over K with PARI
+                # As PARI optimizes the polynomial used to define the number field, we need to precompute
+                # this optimization in Sage and then permute the local invariants correctly
+                x = polygen(QQ, 'x')
+                g = K.pari_polynomial().sage({'x': x})
+                alpha = g.roots(ring=K, multiplicities=False)[0]
+
+                # Compute the number field PARI actually uses, together with isomorphisms to and from K
+                L, to_K, from_K = K.change_generator(alpha)
+
+                # This computation of the permutation relies on the fact that both
+                # Sage and PARI sort real roots of polynomials in increasing order
+                v = K.gen()
+                vals_embed = [sigma(from_K(v)) for sigma in L.embeddings(AA)]
+                perm = Word(vals_embed).standard_permutation()
+
+                # Transfer the primes to PARI and permute the local invariants
+                fin_places_pari = [I.pari_prime() for I in primes]
+                inv_arch_pari = [arg2[i-1] for i in perm]
+
+                # Compute the correct quaternion algebra over L in PARI
+                A = L.__pari__().alginit([2, [fin_places_pari, [QQ(1/2)] * len(fin_places_pari)],
+                                          inv_arch_pari], maxord=0)
+
+                # Obtain representation of A in terms of invariants in L
+                a_L = L(A.algsplittingfield().disc()[1])
+                b_L = L(A.algb())
+
+                # Finally, transfer the result to K
+                a = to_K(a_L)
+                b = to_K(b_L)
+        else:
+            # QuaternionAlgebra(K, a, b)
+            K = arg0
+            a = K(arg1)
+            b = K(arg2)
 
         if not K(2).is_unit():
             raise ValueError("2 is not invertible in %s" % K)