path: root/sagemath-pari-2.13.patch

diff --git a/src/sage/arith/misc.py b/src/sage/arith/misc.py
index 53d9683279..564b99044f 100644
--- a/src/sage/arith/misc.py
+++ b/src/sage/arith/misc.py
@@ -1458,13 +1458,13 @@ def divisors(n):
 
         sage: K.<a> = QuadraticField(7)
         sage: divisors(K.ideal(7))
-        [Fractional ideal (1), Fractional ideal (a), Fractional ideal (7)]
+        [Fractional ideal (1), Fractional ideal (-a), Fractional ideal (7)]
         sage: divisors(K.ideal(3))
         [Fractional ideal (1), Fractional ideal (3),
         Fractional ideal (-a + 2), Fractional ideal (-a - 2)]
         sage: divisors(K.ideal(35))
-        [Fractional ideal (1), Fractional ideal (5), Fractional ideal (a),
-        Fractional ideal (7), Fractional ideal (5*a), Fractional ideal (35)]
+        [Fractional ideal (1), Fractional ideal (5), Fractional ideal (-a),
+        Fractional ideal (7), Fractional ideal (-5*a), Fractional ideal (35)]
 
     TESTS::
 
diff --git a/src/sage/ext_data/pari/simon/ell.gp b/src/sage/ext_data/pari/simon/ell.gp
index 74f0786646..21cff9cbb3 100644
--- a/src/sage/ext_data/pari/simon/ell.gp
+++ b/src/sage/ext_data/pari/simon/ell.gp
@@ -1038,7 +1038,7 @@ if( DEBUGLEVEL_ell >= 1, print(" trivial points on E(K) = ");
   KS2gen = KS2gen[1];
   for( i = 1, #KS2gen,
     KS2gen[i] = nfbasistoalg(bnf, KS2gen[i]));
-  KS2gen = concat(Mod(lift(bnf.tufu),bnf.pol),KS2gen);
+  KS2gen = concat(Mod(lift(concat(bnf.tu[2], bnf.fu)),bnf.pol),KS2gen);
 if( DEBUGLEVEL_ell >= 2,
   print("  #K(b,2)gen          = ",#KS2gen);
   print("  K(b,2)gen = ",KS2gen));
@@ -1072,7 +1072,7 @@ if( DEBUGLEVEL_ell >= 1,
   KS2gen = KS2gen[1];
   for( i = 1, #KS2gen,
     KS2gen[i] = nfbasistoalg(bnf, KS2gen[i]));
-  KS2gen = concat(Mod(lift(bnf.tufu),bnf.pol),KS2gen);
+  KS2gen = concat(Mod(lift(concat(bnf.tu[2], bnf.fu)),bnf.pol),KS2gen);
 if( DEBUGLEVEL_ell >= 2,
   print("  #K(a^2-4b,2)gen     = ",#KS2gen);
   print("  K(a^2-4b,2)gen     = ",KS2gen));
@@ -1244,11 +1244,11 @@ if( DEBUGLEVEL_ell >= 4, print("    bbbnf.clgp = ",bbbnf.clgp));
   SL1 = idealmul(bbbnf,SL0,rnfeltup(rrrnf,bleg));
   SL = idealfactor(bbbnf,SL1)[,1]~;
   sunL = bnfsunit(bbbnf,SL);
-  fondsunL = concat(bbbnf.futu,vector(#sunL[1],i,nfbasistoalg(bbbnf,sunL[1][i])));
+  fondsunL = concat(concat(bbbnf.fu, bbbnf.tu[2]),vector(#sunL[1],i,nfbasistoalg(bbbnf,sunL[1][i])));
   normfondsunL = vector(#fondsunL, i, norm(rnfeltabstorel(rrrnf,fondsunL[i])));
   SK = idealfactor(bnf,idealnorm(bbbnf,SL1))[,1]~;
   sunK = bnfsunit(bnf,SK);
-  fondsunK = concat(bnf.futu,vector(#sunK[1],i,nfbasistoalg(bnf,sunK[1][i])));
+  fondsunK = concat(concat(bnf.fu, bnf.tu[2]),vector(#sunK[1],i,nfbasistoalg(bnf,sunK[1][i])));
   vecbleg = bnfissunit(bnf,sunK,bleg);
   matnorm = matrix(#fondsunK,#normfondsunL,i,j,0);
   for( i = 1, #normfondsunL,
@@ -1345,7 +1345,7 @@ if( DEBUGLEVEL_ell >= 4, print("on factorise bb = ",bb));
       sun = bnfsunit(bnf,idealfactor(bnf,bb)[,1]~);
       fact = lift(bnfissunit(bnf,sun,bb));
 if( DEBUGLEVEL_ell >= 4, print("fact = ",fact));
-      suni = concat(bnf.futu,vector(#sun[1],i,nfbasistoalg(bnf,sun[1][i])));
+      suni = concat(concat(bnf.fu, bnf.tu[2]),vector(#sun[1],i,nfbasistoalg(bnf,sun[1][i])));
       for( i = 1, #suni,
         if( (f = fact[i]>>1), 
           test =0;
@@ -1554,7 +1554,7 @@ if( DEBUGLEVEL_ell >= 3, print("    KS2gen = ",KS2gen[1]));
 
   LS2gen = LS2gen[1];
   LS2 = vector(#LS2gen,i,lift(nfbasistoalg(Lrnf,LS2gen[i])));
-  LS2 = concat(lift(Lrnf.futu),LS2);
+  LS2 = concat(lift(concat(Lrnf.fu, Lrnf.tu[2])),LS2);
 
   LS2 = subst(LS2,'x,ttheta);
   LS2 = LS2*Mod(1,polrel);
@@ -1992,7 +1992,7 @@ if( DEBUGLEVEL_ell >= 2, print("  Algorithm of complete 2-descent"));
   KS2gen = KS2gen[1];
   for( i = 1, #KS2gen,
     KS2gen[i] = nfbasistoalg(bnf, KS2gen[i]));
-  KS2gen = concat(Mod(lift(bnf.tufu),bnf.pol),KS2gen);
+  KS2gen = concat(Mod(lift(concat(bnf.tu[2], bnf.fu)),bnf.pol),KS2gen);
 if( DEBUGLEVEL_ell >= 2,
   print("  #K(S,2)gen = ",#KS2gen);
   print("   K(S,2)gen = ",KS2gen)
diff --git a/src/sage/ext_data/pari/simon/ellQ.gp b/src/sage/ext_data/pari/simon/ellQ.gp
index aede9fc941..27cc124372 100644
--- a/src/sage/ext_data/pari/simon/ellQ.gp
+++ b/src/sage/ext_data/pari/simon/ellQ.gp
@@ -1162,7 +1162,7 @@ if( DEBUGLEVEL_ell >= 4, print("    kerval = ",kerval));
       LS2gen[j]^kerval[j,i]));
 
 \\ Add the units
-  LS2gen = concat(Mod(bnf[8][5],bnf.pol),LS2gen); \\ LS2gen = concat(bnf.fu,LS2gen);
+  LS2gen = concat(bnf.fu,LS2gen); \\ LS2gen = concat(bnf.fu,LS2gen);
 \\ Add also the torsion unit if its order is divisible by p.
   if( bnf[8][4][1]%p == 0, \\ if( bnf.tu[1]%p == 0,
     LS2gen = concat( [Mod(bnf[8][4][2],bnf.pol)], LS2gen)); \\ LS2gen = concat( [Mod(bnf.tu[2],bnf.pol)], LS2gen));
diff --git a/src/sage/geometry/cone.py b/src/sage/geometry/cone.py
index 33b28dbede..1c6873eadd 100644
--- a/src/sage/geometry/cone.py
+++ b/src/sage/geometry/cone.py
@@ -4268,31 +4268,31 @@ class ConvexRationalPolyhedralCone(IntegralRayCollection, Container):
             M(-5,  21,  0, -3),
             M( 0,  -2,  0,  1),
             M(15, -63, 25,  9),
-            M( 2,  -3,  0,  1),
-            M( 1,  -4,  1,  1),
-            M(-1,   3,  0,  0),
             M( 4,  -4,  0,  1),
-            M( 1,  -5,  2,  1),
             M( 3,  -5,  1,  1),
-            M( 6,  -5,  0,  1),
-            M( 3, -13,  5,  2),
             M( 2,  -6,  2,  1),
-            M( 5,  -6,  1,  1),
-            M( 0,   1,  0,  0),
-            M( 8,  -6,  0,  1),
+            M(-1,   7,  0, -1),
+            M( 6, -21,  8,  3),
+            M( 5, -21,  9,  3),
+            M(-1,   3,  0,  0),
+            M( 7, -28, 11,  4),
+            M( 1,  -5,  2,  1),
             M(-2,   8,  0, -1),
+            M( 8,  -6,  0,  1),
+            M( 7,  -7,  1,  1),
+            M( 3, -13,  5,  2),
+            M( 2,  -3,  0,  1),
+            M( 1,  -4,  1,  1),
+            M(-3,  14,  0, -2),
             M(10, -42, 17,  6),
-            M( 7, -28, 11,  4),
-            M( 5, -21,  9,  3),
-            M( 6, -21,  8,  3),
+            M( 1,   0,  0,  0),
+            M( 0,   0,  1,  0),
+            M( 6,  -5,  0,  1),
+            M( 5,  -6,  1,  1),
+            M( 4,  -7,  2,  1),
             M( 5, -14,  5,  2),
             M( 2,  -7,  3,  1),
-            M( 4,  -7,  2,  1),
-            M( 7,  -7,  1,  1),
-            M( 0,   0,  1,  0),
-            M(-3,  14,  0, -2),
-            M(-1,   7,  0, -1),
-            M( 1,   0,  0,  0)
+            M( 0,   1,  0,  0)
             in 4-d lattice M
 
         Not a strictly convex cone::
diff --git a/src/sage/groups/abelian_gps/abelian_group.py b/src/sage/groups/abelian_gps/abelian_group.py
index b61b94c7e4..f7ccab44a3 100644
--- a/src/sage/groups/abelian_gps/abelian_group.py
+++ b/src/sage/groups/abelian_gps/abelian_group.py
@@ -1451,7 +1451,7 @@ class AbelianGroup_class(UniqueRepresentation, AbelianGroupBase):
         EXAMPLES::
 
             sage: AbelianGroup([2,3]).subgroups()
-            [Multiplicative Abelian subgroup isomorphic to C2 x C3 generated by {f0*f1^2},
+            [Multiplicative Abelian subgroup isomorphic to C2 x C3 generated by {f0*f1},
              Multiplicative Abelian subgroup isomorphic to C2 generated by {f0},
              Multiplicative Abelian subgroup isomorphic to C3 generated by {f1},
              Trivial Abelian subgroup]
diff --git a/src/sage/groups/additive_abelian/additive_abelian_group.py b/src/sage/groups/additive_abelian/additive_abelian_group.py
index 3dcedeb7e3..48747e4870 100644
--- a/src/sage/groups/additive_abelian/additive_abelian_group.py
+++ b/src/sage/groups/additive_abelian/additive_abelian_group.py
@@ -61,18 +61,18 @@ def AdditiveAbelianGroup(invs, remember_generators = True):
         ((1, 0, 0), (0, 1, 0), (0, 0, 1))
         sage: [H.0, H.1, H.2]
         [(1, 0, 0), (0, 1, 0), (0, 0, 1)]
-        sage: p=H.0+H.1+6*H.2; p
-        (1, 1, 6)
+        sage: p=2*H.0+H.1+6*H.2; p
+        (2, 1, 6)
 
         sage: H.smith_form_gens()
-        ((2, 1, 0), (0, 0, 1))
+        ((1, 1, 0), (0, 0, 1))
         sage: q=H.linear_combination_of_smith_form_gens([5,6]); q
-        (1, 1, 6)
+        (2, 1, 6)
         sage: p==q
         True
 
-        sage: r=H(vector([1,1,6])); r
-        (1, 1, 6)
+        sage: r=H(vector([2,1,6])); r
+        (2, 1, 6)
         sage: p==r
         True
 
@@ -85,21 +85,21 @@ def AdditiveAbelianGroup(invs, remember_generators = True):
 
         sage: G=AdditiveAbelianGroup([3,2,0], remember_generators=False)
         sage: G.gens()
-        ((2, 1, 0), (0, 0, 1))
+        ((1, 1, 0), (0, 0, 1))
         sage: [G.0, G.1]
-        [(2, 1, 0), (0, 0, 1)]
+        [(1, 1, 0), (0, 0, 1)]
         sage: p=5*G.0+6*G.1; p
-        (1, 1, 6)
+        (2, 1, 6)
 
         sage: H.smith_form_gens()
-        ((2, 1, 0), (0, 0, 1))
+        ((1, 1, 0), (0, 0, 1))
         sage: q=G.linear_combination_of_smith_form_gens([5,6]); q
-        (1, 1, 6)
+        (2, 1, 6)
         sage: p==q
         True
 
-        sage: r=G(vector([1,1,6])); r
-        (1, 1, 6)
+        sage: r=G(vector([2,1,6])); r
+        (2, 1, 6)
         sage: p==r
         True
 
@@ -427,7 +427,7 @@ class AdditiveAbelianGroup_fixed_gens(AdditiveAbelianGroup_class):
             sage: G.gens()
             ((1, 0), (0, 1))
             sage: G.smith_form_gens()
-            ((1, 2),)
+            ((1, 1),)
         """
         return self._orig_gens
 
diff --git a/src/sage/groups/fqf_orthogonal.py b/src/sage/groups/fqf_orthogonal.py
index 3c5190589b..c536f03064 100644
--- a/src/sage/groups/fqf_orthogonal.py
+++ b/src/sage/groups/fqf_orthogonal.py
@@ -156,9 +156,9 @@ class FqfOrthogonalGroup(AbelianGroupAutomorphismGroup_subgroup):
         sage: S2 = 9 * T
         sage: Q = S1/S2
         sage: G = T.orthogonal_group()
-        sage: g = G(matrix(ZZ, 2, [8, 0, 0, 1]))
+        sage: g = G(matrix(ZZ, 2, [7, 12, 8, 19]))
         sage: Q.1 * g
-        (0, 2)
+        (0, 1)
     """
     Element = FqfIsometry
 
diff --git a/src/sage/lfunctions/pari.py b/src/sage/lfunctions/pari.py
index f810157b3e..f621b57c67 100644
--- a/src/sage/lfunctions/pari.py
+++ b/src/sage/lfunctions/pari.py
@@ -421,7 +421,7 @@ class LFunction(SageObject):
         sage: L.derivative(1,E.rank())
         1.51863300057685
         sage: L.taylor_series(1,4)
-        -3...e-19 + (...e-19)*z + 0.759316500288427*z^2 - 0.430302337583362*z^3 + O(z^4)
+        ...e-19 + (...e-19)*z + 0.759316500288427*z^2 - 0.430302337583362*z^3 + O(z^4)
 
     .. RUBRIC:: Number field
 
diff --git a/src/sage/libs/pari/__init__.py b/src/sage/libs/pari/__init__.py
index 77eda66097..3fa4618631 100644
--- a/src/sage/libs/pari/__init__.py
+++ b/src/sage/libs/pari/__init__.py
@@ -161,12 +161,12 @@ exact object. Therefore, you should set the precision for each method
 call individually::
 
     sage: e = pari([0,0,0,-82,0]).ellinit()
-    sage: eta1 = e.elleta(precision=100)[0]
+    sage: eta1 = e.elleta(precision=50)[0]
     sage: eta1.sage()
     3.6054636014326520859158205642077267748
     sage: eta1 = e.elleta(precision=180)[0]
     sage: eta1.sage()
-    3.60546360143265208591582056420772677481026899659802474544
+    3.605463601432652085915820564207726774810268996598024745444380641429820491740
 
 """
 
diff --git a/src/sage/libs/pari/convert_sage.pyx b/src/sage/libs/pari/convert_sage.pyx
index 33aa4d2d46..b506227cff 100644
--- a/src/sage/libs/pari/convert_sage.pyx
+++ b/src/sage/libs/pari/convert_sage.pyx
@@ -144,6 +144,20 @@ cpdef gen_to_sage(Gen z, locals=None):
         sage: a.parent()
         Complex Field with 64 bits of precision
 
+        sage: z = pari('1 + 1.0*I'); z
+        1 + 1.00000000000000*I
+        sage: a = gen_to_sage(z); a
+        1.00000000000000000 + 1.00000000000000000*I
+        sage: a.parent()
+        Complex Field with 64 bits of precision
+
+        sage: z = pari('1.0 + 1*I'); z
+        1.00000000000000 + I
+        sage: a = gen_to_sage(z); a
+        1.00000000000000000 + 1.00000000000000000*I
+        sage: a.parent()
+        Complex Field with 64 bits of precision
+
     Converting polynomials::
 
         sage: f = pari('(2/3)*x^3 + x - 5/7 + y')
@@ -241,7 +255,9 @@ cpdef gen_to_sage(Gen z, locals=None):
     elif t == t_FRAC:
         return Rational(z)
     elif t == t_REAL:
-        prec = prec_words_to_bits(z.precision())
+        prec = z.bitprecision()
+        if prec.type() == 't_INFINITY':
+            prec = 53
         return RealField(prec)(z)
     elif t == t_COMPLEX:
         real = z.real()
@@ -251,14 +267,19 @@ cpdef gen_to_sage(Gen z, locals=None):
         if tx in [t_INTMOD, t_PADIC] or ty in [t_INTMOD, t_PADIC]:
             raise NotImplementedError("No conversion to python available for t_COMPLEX with t_INTMOD or t_PADIC components")
         if tx == t_REAL or ty == t_REAL:
-            xprec = real.precision()  # will be 0 if exact
-            yprec = imag.precision()  # will be 0 if exact
-            if xprec == 0:
-                prec = prec_words_to_bits(yprec)
-            elif yprec == 0:
-                prec = prec_words_to_bits(xprec)
+            xprec = real.bitprecision()  # will be 0 if exact
+            yprec = imag.bitprecision()  # will be 0 if exact
+            if xprec == 0 or yprec == 0:
+                raise RuntimeError
+            if xprec.type() == 't_INFINITY':
+                if yprec.type() == 't_INFINITY':
+                    prec = 53
+                else:
+                    prec = yprec
+            elif yprec.type() == 't_INFINITY':
+                prec = xprec
             else:
-                prec = max(prec_words_to_bits(xprec), prec_words_to_bits(yprec))
+                prec = max(xprec, yprec)
 
             R = RealField(prec)
             C = ComplexField(prec)
diff --git a/src/sage/libs/pari/tests.py b/src/sage/libs/pari/tests.py
index eb1c04f1e6..40329aea82 100644
--- a/src/sage/libs/pari/tests.py
+++ b/src/sage/libs/pari/tests.py
@@ -135,7 +135,7 @@ Some more exotic examples::
 
     sage: K.<a> = NumberField(polygen(QQ)^3 - 2)
     sage: pari(K)
-    [y^3 - 2, [1, 1], -108, 1, [[1, 1.25992104989487, 1.58740105196820; 1, -0.629960524947437 + 1.09112363597172*I, -0.793700525984100 - 1.37472963699860*I], [1, 1.25992104989487, 1.58740105196820; 1, 0.461163111024285, -2.16843016298270; 1, -1.72108416091916, 0.581029111014503], [1, 1, 2; 1, 0, -2; 1, -2, 1], [3, 0, 0; 0, 0, 6; 0, 6, 0], [6, 0, 0; 0, 6, 0; 0, 0, 3], [2, 0, 0; 0, 0, 1; 0, 1, 0], [2, [0, 0, 2; 1, 0, 0; 0, 1, 0]], []], [1.25992104989487, -0.629960524947437 + 1.09112363597172*I], [1, y, y^2], [1, 0, 0; 0, 1, 0; 0, 0, 1], [1, 0, 0, 0, 0, 2, 0, 2, 0; 0, 1, 0, 1, 0, 0, 0, 0, 2; 0, 0, 1, 0, 1, 0, 1, 0, 0]]
+    [y^3 - 2, [1, 1], -108, 1, [[1, 1.25992104989487, 1.58740105196820; 1, -0.629960524947437 + 1.09112363597172*I, -0.793700525984100 - 1.37472963699860*I], [1, 1.25992104989487, 1.58740105196820; 1, 0.461163111024285, -2.16843016298270; 1, -1.72108416091916, 0.581029111014503], [16, 20, 25; 16, 7, -35; 16, -28, 9], [3, 0, 0; 0, 0, 6; 0, 6, 0], [6, 0, 0; 0, 6, 0; 0, 0, 3], [2, 0, 0; 0, 0, 1; 0, 1, 0], [2, [0, 0, 2; 1, 0, 0; 0, 1, 0]], [2, 3]], [1.25992104989487, -0.629960524947437 + 1.09112363597172*I], [1, y, y^2], [1, 0, 0; 0, 1, 0; 0, 0, 1], [1, 0, 0, 0, 0, 2, 0, 2, 0; 0, 1, 0, 1, 0, 0, 0, 0, 2; 0, 0, 1, 0, 1, 0, 1, 0, 0]]
 
     sage: E = EllipticCurve('37a1')
     sage: pari(E)
@@ -375,13 +375,13 @@ Constructors::
     sage: pari('["bc","ab","bc"]').Set()
     ["ab", "bc"]
 
-    sage: pari([65,66,123]).Strchr()
+    sage: pari([65,66,123]).strchr()
     "AB{"
     sage: pari('"Sage"').Vecsmall()
     Vecsmall([83, 97, 103, 101])
-    sage: _.Strchr()
+    sage: _.strchr()
     "Sage"
-    sage: pari([83, 97, 103, 101]).Strchr()
+    sage: pari([83, 97, 103, 101]).strchr()
     "Sage"
 
 Basic functions::
@@ -448,7 +448,7 @@ Basic functions::
     sage: pari('x').component(0)
     Traceback (most recent call last):
     ...
-    PariError: non-existent component: index < 1
+    PariError: nonexistent component: index < 1
 
     sage: pari('x+1').conj()
     x + 1
@@ -767,7 +767,7 @@ Transcendental functions::
     sage: pari(2).besseli(3+i)
     1.12539407613913 + 2.08313822670661*I
     sage: C.<i> = ComplexField()
-    sage: pari(2+i).besseln(3)
+    sage: pari(2+i).bessely(3)
     -0.280775566958244 - 0.486708533223726*I
 
     sage: pari(1.5).cos()
@@ -822,7 +822,7 @@ Transcendental functions::
     sage: pari(-1).gamma()
     Traceback (most recent call last):
     ...
-    PariError: domain error in gamma: argument = non-positive integer
+    PariError: domain error in gamma: argument = nonpositive integer
 
     sage: pari(2).gammah()
     1.32934038817914
@@ -1633,7 +1633,7 @@ General number fields::
 
     sage: x = QQ['x'].0; nf = pari(x^2 + 2).nfinit()
     sage: nf.nfgaloisconj()
-    [x, -x]~
+    [-x, x]~
     sage: nf = pari(x^3 + 2).nfinit()
     sage: nf.nfgaloisconj()
     [x]~
@@ -1676,7 +1676,7 @@ General number fields::
     [[1, [7605, 4]~, [5610, 5]~, [7913, -6]~; 0, 1, 0, -1; 0, 0, 1, 0; 0, 0, 0, 1], [[19320, 13720; 0, 56], [2, 1; 0, 1], 1, 1]]
 
     sage: pari('x^3 - 17').nfinit()
-    [x^3 - 17, [1, 1], -867, 3, [[1, 1.68006914259990, 2.57128159065824; 1, -0.340034571299952 - 2.65083754153991*I, -1.28564079532912 + 2.22679517779329*I], [1, 1.68006914259990, 2.57128159065824; 1, -2.99087211283986, 0.941154382464174; 1, 2.31080297023995, -3.51243597312241], [1, 2, 3; 1, -3, 1; 1, 2, -4], [3, 1, 0; 1, -11, 17; 0, 17, 0], [51, 0, 16; 0, 17, 3; 0, 0, 1], [17, 0, -1; 0, 0, 3; -1, 3, 2], [51, [-17, 6, -1; 0, -18, 3; 1, 0, -16]], [3, 17]], [2.57128159065824, -1.28564079532912 + 2.22679517779329*I], [3, x^2 - x + 1, 3*x], [1, 0, -1; 0, 0, 3; 0, 1, 1], [1, 0, 0, 0, -4, 6, 0, 6, -1; 0, 1, 0, 1, 1, -1, 0, -1, 3; 0, 0, 1, 0, 2, 0, 1, 0, 1]]
+    [x^3 - 17, [1, 1], -867, 3, [[1, 1.68006914259990, 2.57128159065824; 1, -0.340034571299952 - 2.65083754153991*I, -1.28564079532912 + 2.22679517779329*I], [1, 1.68006914259990, 2.57128159065824; 1, -2.99087211283986, 0.941154382464174; 1, 2.31080297023995, -3.51243597312241], [16, 27, 41; 16, -48, 15; 16, 37, -56], [3, 1, 0; 1, -11, 17; 0, 17, 0], [51, 0, 16; 0, 17, 3; 0, 0, 1], [17, 0, -1; 0, 0, 3; -1, 3, 2], [51, [-17, 6, -1; 0, -18, 3; 1, 0, -16]], [3, 17]], [2.57128159065824, -1.28564079532912 + 2.22679517779329*I], [3, x^2 - x + 1, 3*x], [1, 0, -1; 0, 0, 3; 0, 1, 1], [1, 0, 0, 0, -4, 6, 0, 6, -1; 0, 1, 0, 1, 1, -1, 0, -1, 3; 0, 0, 1, 0, 2, 0, 1, 0, 1]]
     sage: pari('x^2 + 10^100 + 1').nfinit()
     [...]
     sage: pari('1.0').nfinit()
@@ -1737,7 +1737,7 @@ General number fields::
     sage: pari(-23).quadhilbert()
     x^3 - x^2 + 1
     sage: pari(145).quadhilbert()
-    x^4 - 6*x^2 - 5*x - 1
+    x^4 - x^3 - 5*x^2 - x + 1
     sage: pari(-12).quadhilbert()   # Not fundamental
     Traceback (most recent call last):
     ...
@@ -1762,7 +1762,7 @@ library::
     sage: e = pari([0,0,0,-82,0]).ellinit()
     sage: eta1 = e.elleta(precision=100)[0]
     sage: eta1.sage()
-    3.6054636014326520859158205642077267748
+    3.60546360143265208591582056420772677481026899659802474544
     sage: eta1 = e.elleta(precision=180)[0]
     sage: eta1.sage()
     3.60546360143265208591582056420772677481026899659802474544
diff --git a/src/sage/matrix/matrix1.pyx b/src/sage/matrix/matrix1.pyx
index 0455c169fe..4fa3cb9b2d 100644
--- a/src/sage/matrix/matrix1.pyx
+++ b/src/sage/matrix/matrix1.pyx
@@ -87,7 +87,7 @@ cdef class Matrix(Matrix0):
             [1.000000000, 2.000000000; 3.000000000, 1.000000000] # 32-bit
             [1.00000000000000, 2.00000000000000; 3.00000000000000, 1.00000000000000] # 64-bit
             sage: b[0][0].precision()    # in words
-            3
+            19
         """
         from sage.libs.pari.all import pari
         return pari.matrix(self._nrows, self._ncols, self._list())
diff --git a/src/sage/modular/local_comp/liftings.py b/src/sage/modular/local_comp/liftings.py
index 91bfd2448d..f9bf800b22 100644
--- a/src/sage/modular/local_comp/liftings.py
+++ b/src/sage/modular/local_comp/liftings.py
@@ -222,9 +222,9 @@ def lift_for_SL(A, N=None):
     TESTS::
 
         sage: lift_for_SL(matrix(3,3,[1,2,0,3,4,0,0,0,1]),3)
-        [10 14  3]
-        [ 9 10  3]
-        [ 3  3  1]
+        [ -2  -4   3]
+        [ -9 -14   3]
+        [ -6  -9   1]
 
         sage: A = matrix(Zmod(7), 2, [1,0,0,1])
         sage: L = lift_for_SL(A)
diff --git a/src/sage/modular/local_comp/smoothchar.py b/src/sage/modular/local_comp/smoothchar.py
index 6dedb26e7e..3b1ac90460 100644
--- a/src/sage/modular/local_comp/smoothchar.py
+++ b/src/sage/modular/local_comp/smoothchar.py
@@ -1617,8 +1617,8 @@ class SmoothCharacterGroupRamifiedQuadratic(SmoothCharacterGroupGeneric):
             sage: G = SmoothCharacterGroupRamifiedQuadratic(3, 1, QQ)
             sage: s = G.number_field().gen()
             sage: G.discrete_log(4, 3 + 2*s)
-            [5, 1, 1, 1]
-            sage: gs = G.unit_gens(4); gs[0]^5 * gs[1] * gs[2] * gs[3] - (3 + 2*s) in G.ideal(4)
+            [1, 2, 2, 1]
+            sage: gs = G.unit_gens(4); gs[0] * gs[1]^2 * gs[2]^2 * gs[3] - (3 + 2*s) in G.ideal(4)
             True
         """
         x = self.number_field().coerce(x)
diff --git a/src/sage/modular/modsym/p1list_nf.py b/src/sage/modular/modsym/p1list_nf.py
index 7a91d353e3..b71fc8aac4 100644
--- a/src/sage/modular/modsym/p1list_nf.py
+++ b/src/sage/modular/modsym/p1list_nf.py
@@ -956,7 +956,7 @@ class P1NFList(SageObject):
             sage: N = k.ideal(a + 1)
             sage: P = P1NFList(N)
             sage: u = k.unit_group().gens_values(); u
-            [-1, a^3 + a^2 + a + 12, a^3 + 3*a^2 - 1]
+            [-1, -a^3 - a^2 - a - 12, -a^3 - 3*a^2 + 1]
             sage: P.apply_J_epsilon(3, u[2]^2)==P.apply_J_epsilon(P.apply_J_epsilon(3, u[2]),u[2])
             True
         """
diff --git a/src/sage/modular/multiple_zeta.py b/src/sage/modular/multiple_zeta.py
index d42be82bd9..5f0057474c 100644
--- a/src/sage/modular/multiple_zeta.py
+++ b/src/sage/modular/multiple_zeta.py
@@ -457,7 +457,7 @@ class MultizetaValues(UniqueRepresentation):
         """
         self.prec = int(prec)
         self.max_weight = int(max_weight)
-        self._data = pari.zetamultall(self.max_weight, self.prec)
+        self._data = pari.zetamultall(self.max_weight, precision=self.prec)
 
     def update(self, max_weight, prec):
         """
diff --git a/src/sage/modules/fg_pid/fgp_element.py b/src/sage/modules/fg_pid/fgp_element.py
index 53857d31a0..6557b7f308 100644
--- a/src/sage/modules/fg_pid/fgp_element.py
+++ b/src/sage/modules/fg_pid/fgp_element.py
@@ -39,7 +39,7 @@ class FGP_Element(ModuleElement):
         sage: V = span([[1/2,1,1],[3/2,2,1],[0,0,1]],ZZ); W = V.span([2*V.0+4*V.1, 9*V.0+12*V.1, 4*V.2])
         sage: Q = V/W
         sage: x = Q(V.0-V.1); x #indirect doctest
-        (0, 3)
+        (0, 9)
         sage: isinstance(x, sage.modules.fg_pid.fgp_element.FGP_Element)
         True
         sage: type(x)
@@ -94,14 +94,14 @@ class FGP_Element(ModuleElement):
             sage: Q.1
             (0, 1)
             sage: Q.0.lift()
-            (0, 0, 1)
+            (0, -6, 1)
             sage: Q.1.lift()
-            (0, 2, 0)
+            (0, -2, 0)
             sage: x = Q(V.0); x
-            (0, 4)
+            (0, 8)
             sage: x.lift()
             (1/2, 0, 0)
-            sage: x == 4*Q.1
+            sage: x == 8*Q.1
             True
             sage: x.lift().parent() == V
             True
@@ -158,9 +158,9 @@ class FGP_Element(ModuleElement):
         We test canonical coercion from V and W.
 
             sage: Q.0 + V.0
-            (1, 4)
+            (1, 8)
             sage: V.0 + Q.0
-            (1, 4)
+            (1, 8)
             sage: W.0 + Q.0
             (1, 0)
             sage: W.0 + Q.0 == Q.0
@@ -291,7 +291,7 @@ class FGP_Element(ModuleElement):
             sage: V = span([[1/2,1,1],[3/2,2,1],[0,0,1]],ZZ); W = V.span([2*V.0+4*V.1, 9*V.0+12*V.1, 4*V.2])
             sage: Q = V/W
             sage: Q(V.1)._repr_()
-            '(0, 1)'
+            '(0, 11)'
         """
         return repr(self.vector())
 
diff --git a/src/sage/modules/fg_pid/fgp_module.py b/src/sage/modules/fg_pid/fgp_module.py
index 35f7cdaa17..5a700d808e 100644
--- a/src/sage/modules/fg_pid/fgp_module.py
+++ b/src/sage/modules/fg_pid/fgp_module.py
@@ -70,17 +70,17 @@ the technical note has a V that need not be equal to V0, in general. ::
     sage: M0.optimized()[0].V()
     Free module of degree 3 and rank 2 over Integer Ring
     User basis matrix:
-    [0 0 1]
-    [0 2 0]
+    [ 0 -8  1]
+    [ 0 -2  0]
 
 Create elements of M0 either by coercing in elements of V0, getting generators,
 or coercing in a list or tuple or coercing in 0. Finally, one can express an
 element as a linear combination of the smith form generators ::
 
     sage: M0(V0.0)
-    (0, 14)
+    (0, 2)
     sage: M0(V0.0 + W0.0)  # no difference modulo W0
-    (0, 14)
+    (0, 2)
     sage: M0.linear_combination_of_smith_form_gens([3,20])
     (3, 4)
     sage: 3*M0.0 + 20*M0.1
@@ -93,9 +93,9 @@ coerces to V0, then take the equivalence class modulo W0. ::
     sage: x = M0.0 - M0.1; x
     (1, 15)
     sage: x.lift()
-    (0, -2, 1)
+    (0, -6, 1)
     sage: M0(vector([1/2,0,0]))
-    (0, 14)
+    (0, 2)
     sage: x.additive_order()
     16
 
@@ -144,7 +144,7 @@ You can explicitly coerce elements of the kernel into M0 though. ::
     sage: M0(K.0)
     (2, 0)
     sage: M0(K.1)
-    (3, 1)
+    (1, 13)
     sage: f(M0(K.0))
     (0)
     sage: f(M0(K.1))
@@ -180,7 +180,7 @@ TESTS::
     sage: Q.linear_combination_of_smith_form_gens([1,3])
     (1, 3)
     sage: Q(V([1,3,4]))
-    (0, 11)
+    (0, 1)
     sage: Q(W([1,16,0]))
     (0, 0)
     sage: V = span([[1/2,1,1],[3/2,2,1],[0,0,1]],QQ)
@@ -632,7 +632,7 @@ class FGP_Module_class(Module):
             sage: W = V.span([2*V.0+4*V.1, 9*V.0+12*V.1, 4*V.2])
             sage: Q = V/W
             sage: x = Q(V.0-V.1); x  # indirect doctest
-            (0, 3)
+            (0, 9)
             sage: type(x)
             <class 'sage.modules.fg_pid.fgp_module.FGP_Module_class_with_category.element_class'>
             sage: x is Q(x)
@@ -931,9 +931,9 @@ class FGP_Module_class(Module):
             sage: Q = V/W
             sage: Q._smith_form()
             (
-            [ 1  0  0]  [1 0 0]  [ 1  0 -8]
-            [ 0  4  0]  [0 0 1]  [ 0  0  1]
-            [ 0  0 12], [0 1 0], [ 0  1  0]
+            [ 1  0  0]  [ 1  0  1]  [ 1 -4 20]
+            [ 0  4  0]  [ 0 -1  1]  [ 0  0 -1]
+            [ 0  0 12], [ 0 -1  0], [ 0  1 -3]
             )
         """
         return self._relative_matrix().smith_form()
@@ -1026,7 +1026,7 @@ class FGP_Module_class(Module):
             sage: Q.smith_form_gens()
             ((1, 0), (0, 1))
             sage: [x.lift() for x in Q.smith_form_gens()]
-            [(0, 0, 1), (0, 1, 0)]
+            [(0, -3, 1), (0, -1, 0)]
         """
         # Get the rightmost transformation in the Smith form
         _, _, X = self._smith_form()
@@ -1067,18 +1067,18 @@ class FGP_Module_class(Module):
             [  0   0   0 1/3   0]
             [  0   0   0   0 2/3]
             sage: D.gens_to_smith()
-            [0 3 0]
+            [0 3 6]
             [0 0 3]
-            [0 2 0]
-            [1 0 0]
-            [0 0 4]
+            [0 4 4]
+            [1 2 0]
+            [0 0 8]
             sage: T = D.gens_to_smith()*D.smith_to_gens()
             sage: T
-            [ 3  0 15  0  0]
-            [ 0 33  0  0  3]
-            [ 2  0 10  0  0]
-            [ 0  0  0  1  0]
-            [ 0 44  0  0  4]
+            [27 48  3  0 60]
+            [12 21  0  0 24]
+            [20 36  4  0 48]
+            [ 2  4  3  1  9]
+            [32 56  0  0 64]
 
         The matrix `T` now satisfies a certain congruence::
 
@@ -1120,14 +1120,14 @@ class FGP_Module_class(Module):
             [  0   0   0 1/3   0]
             [  0   0   0   0 2/3]
             sage: D.smith_to_gens()
-            [ 0  0  0  1  0]
-            [ 1  0  5  0  0]
-            [ 0 11  0  0  1]
+            [0 0 1 1 1]
+            [1 2 1 0 4]
+            [4 7 0 0 8]
             sage: T = D.smith_to_gens()*D.gens_to_smith()
             sage: T
-            [ 1  0  0]
-            [ 0 13  0]
-            [ 0  0 37]
+            [  1   6  12]
+            [  0   7  48]
+            [  0  12 109]
 
         This matrix satisfies the congruence::
 
@@ -1148,7 +1148,7 @@ class FGP_Module_class(Module):
         of the user defined generators that is x::
 
             sage: x.vector() * D.smith_to_gens()
-            (2, 33, 10, 1, 3)
+            (14, 25, 3, 1, 33)
         """
         if self.base_ring() != ZZ:
             # it is not
@@ -1196,7 +1196,7 @@ class FGP_Module_class(Module):
              sage: gens = [V(g) for g in gens]
              sage: D = FGP_with_gens(V, W, gens)
              sage: D.gens()
-             ((0, 3, 0), (0, 0, 3), (0, 2, 0), (1, 0, 0), (0, 0, 8))
+             ((0, 3, 6), (0, 0, 9), (0, 4, 4), (1, 2, 0), (0, 0, 4))
 
 
         We create some element of D::
@@ -1209,12 +1209,12 @@ class FGP_Module_class(Module):
 
             sage: v = D.gens_vector(x)
             sage: v
-            (2, 9, 10, 1, 33)
+            (26, 11, 3, 1, 18)
 
         The output can be further reduced::
 
             sage: D.gens_vector(x, reduce=True)
-            (0, 1, 1, 1, 0)
+            (0, 3, 0, 1, 0)
 
         Let us check::
 
@@ -1278,28 +1278,28 @@ class FGP_Module_class(Module):
             sage: O.V()
             Free module of degree 3 and rank 2 over Integer Ring
             User basis matrix:
-            [0 0 1]
-            [0 2 0]
+            [ 0 -6  1]
+            [ 0 -2  0]
             sage: phi = Q.hom([Q.0, 4*Q.1])
             sage: x = Q(V.0); x
-            (0, 4)
+            (0, 8)
             sage: Q.coordinate_vector(x, reduce=True)
-            (0, 4)
+            (0, 8)
             sage: Q.coordinate_vector(-x, reduce=False) # random
-            (0, -4)
-            sage: x == 4*Q.1
+            (0, -8)
+            sage: x == 8*Q.1
             True
             sage: x = Q(V.1); x
-            (0, 1)
+            (0, 11)
             sage: Q.coordinate_vector(x)
-            (0, 1)
-            sage: x == Q.1
+            (0, -1)
+            sage: x == -Q.1
             True
             sage: x = Q(V.2); x
-            (1, 0)
+            (1, 9)
             sage: Q.coordinate_vector(x)
-            (1, 0)
-            sage: x == Q.0
+            (1, -3)
+            sage: x == Q.0-3*Q.1
             True
         """
         try:
@@ -1407,8 +1407,8 @@ class FGP_Module_class(Module):
             sage: O.V()
             Free module of degree 3 and rank 2 over Integer Ring
             User basis matrix:
-            [0 0 1]
-            [0 1 0]
+            [ 0 -3  1]
+            [ 0 -1  0]
             sage: O.W()
             Free module of degree 3 and rank 2 over Integer Ring
             Echelon basis matrix:
@@ -1699,7 +1699,7 @@ class FGP_Module_class(Module):
             sage: V = span([[1/2,1,1],[3/2,2,1],[0,0,1]],ZZ); W = V.span([2*V.0+4*V.1, 9*V.0+12*V.1, 4*V.2])
             sage: Q = V/W
             sage: Q.random_element()
-            (1, 10)
+            (1, 11)
         """
         return self(self._V.random_element(*args, **kwds))
 
diff --git a/src/sage/modules/fg_pid/fgp_morphism.py b/src/sage/modules/fg_pid/fgp_morphism.py
index fbbd479c6e..679965df96 100644
--- a/src/sage/modules/fg_pid/fgp_morphism.py
+++ b/src/sage/modules/fg_pid/fgp_morphism.py
@@ -258,20 +258,20 @@ class FGP_Morphism(Morphism):
             sage: O.V()
             Free module of degree 3 and rank 2 over Integer Ring
             User basis matrix:
-            [0 0 1]
-            [0 2 0]
+            [ 0 -6  1]
+            [ 0 -2  0]
             sage: phi = Q.hom([Q.0, 4*Q.1])
             sage: x = Q(V.0); x
-            (0, 4)
-            sage: x == 4*Q.1
+            (0, 8)
+            sage: x == 8*Q.1
             True
             sage: x in O.V()
             False
             sage: phi(x)
-            (0, 4)
+            (0, 8)
             sage: phi(4*Q.1)
             (0, 4)
-            sage: phi(4*Q.1) == phi(x)
+            sage: phi(8*Q.1) == phi(x)
             True
         """
         from .fgp_module import is_FGP_Module
diff --git a/src/sage/modules/torsion_quadratic_module.py b/src/sage/modules/torsion_quadratic_module.py
index 39e7065ac4..6086d1180d 100644
--- a/src/sage/modules/torsion_quadratic_module.py
+++ b/src/sage/modules/torsion_quadratic_module.py
@@ -1230,24 +1230,24 @@ class TorsionQuadraticModule(FGP_Module_class, CachedRepresentation):
             sage: q.twist(-1)
             Finite quadratic module over Integer Ring with invariants (3, 9)
             Gram matrix of the quadratic form with values in Q/Z:
-            [2/3   0]
-            [  0 8/9]
+            [2/3 1/3]
+            [1/3 8/9]
 
         This form is defined modulo `3`::
 
             sage: q.twist(3)
             Finite quadratic module over Integer Ring with invariants (3, 9)
             Gram matrix of the quadratic form with values in Q/3Z:
-            [  1   0]
-            [  0 1/3]
+            [  1   2]
+            [  2 1/3]
 
         The next form is defined modulo `4`::
 
             sage: q.twist(4)
             Finite quadratic module over Integer Ring with invariants (3, 9)
             Gram matrix of the quadratic form with values in Q/4Z:
-            [4/3   0]
-            [  0 4/9]
+            [4/3 8/3]
+            [8/3 4/9]
         """
         s = self.base_ring().fraction_field()(s)
         n = self.V().degree()
diff --git a/src/sage/quadratic_forms/genera/genus.py b/src/sage/quadratic_forms/genera/genus.py
index 2e773616e2..0daa55f84f 100644
--- a/src/sage/quadratic_forms/genera/genus.py
+++ b/src/sage/quadratic_forms/genera/genus.py
@@ -2897,18 +2897,18 @@ class GenusSymbol_global_ring(object):
             sage: GS.discriminant_form()
             Finite quadratic module over Integer Ring with invariants (2, 2, 4, 24)
             Gram matrix of the quadratic form with values in Q/2Z:
-            [ 1/2    0    0    0]
-            [   0  3/2    0    0]
-            [   0    0  7/4    0]
-            [   0    0    0 7/24]
+            [  1/2     0     0     0]
+            [    0   1/2   1/2     0]
+            [    0   1/2   7/4     0]
+            [    0     0     0 31/24]
             sage: A = matrix.diagonal(ZZ, [1, -4, 6, 8])
             sage: GS = Genus(A)
             sage: GS.discriminant_form()
             Finite quadratic module over Integer Ring with invariants (2, 4, 24)
             Gram matrix of the quadratic form with values in Q/Z:
-            [ 1/2    0    0]
-            [   0  3/4    0]
-            [   0    0 7/24]
+            [  1/2   1/2     0]
+            [  1/2   3/4     0]
+            [    0     0 19/24]
         """
         from sage.modules.torsion_quadratic_module import TorsionQuadraticForm
         qL = []
diff --git a/src/sage/rings/finite_rings/finite_field_constructor.py b/src/sage/rings/finite_rings/finite_field_constructor.py
index 81ff16b1b6..46d56b923f 100644
--- a/src/sage/rings/finite_rings/finite_field_constructor.py
+++ b/src/sage/rings/finite_rings/finite_field_constructor.py
@@ -284,11 +284,6 @@ class FiniteFieldFactory(UniqueFactory):
     (a generator of the multiplicative group), use
     ``modulus="primitive"`` if you need this::
 
-        sage: K.<a> = GF(5^40)
-        sage: a.multiplicative_order()
-        189478062869360049565633138
-        sage: a.is_square()
-        True
         sage: K.<b> = GF(5^40, modulus="primitive")
         sage: b.multiplicative_order()
         9094947017729282379150390624
diff --git a/src/sage/rings/finite_rings/integer_mod_ring.py b/src/sage/rings/finite_rings/integer_mod_ring.py
index 0dcef0d21a..7b4f8b1b4b 100644
--- a/src/sage/rings/finite_rings/integer_mod_ring.py
+++ b/src/sage/rings/finite_rings/integer_mod_ring.py
@@ -626,7 +626,7 @@ class IntegerModRing_generic(quotient_ring.QuotientRing_generic):
             sage: Integers(5).multiplicative_subgroups()
             ((2,), (4,), ())
             sage: Integers(15).multiplicative_subgroups()
-            ((11, 7), (4, 11), (8,), (11,), (14,), (7,), (4,), ())
+            ((14, 13), (4, 11), (8,), (11,), (14,), (7,), (4,), ())
             sage: Integers(2).multiplicative_subgroups()
             ((),)
             sage: len(Integers(341).multiplicative_subgroups())
diff --git a/src/sage/rings/finite_rings/residue_field.pyx b/src/sage/rings/finite_rings/residue_field.pyx
index cd4c2212c3..007a68e4c0 100644
--- a/src/sage/rings/finite_rings/residue_field.pyx
+++ b/src/sage/rings/finite_rings/residue_field.pyx
@@ -1055,7 +1055,7 @@ cdef class ReductionMap(Map):
             sage: f = k.convert_map_from(K)
             sage: s = f.section(); s
             Lifting map:
-              From: Residue field in abar of Fractional ideal (14*a^4 - 24*a^3 - 26*a^2 + 58*a - 15)
+              From: Residue field in abar of Fractional ideal (-14*a^4 + 24*a^3 + 26*a^2 - 58*a + 15)
               To:   Number Field in a with defining polynomial x^5 - 5*x + 2
             sage: s(k.gen())
             a
@@ -1268,7 +1268,7 @@ cdef class ResidueFieldHomomorphism_global(RingHomomorphism):
             sage: f = k.coerce_map_from(K.ring_of_integers())
             sage: s = f.section(); s
             Lifting map:
-              From: Residue field in abar of Fractional ideal (14*a^4 - 24*a^3 - 26*a^2 + 58*a - 15)
+              From: Residue field in abar of Fractional ideal (-14*a^4 + 24*a^3 + 26*a^2 - 58*a + 15)
               To:   Maximal Order in Number Field in a with defining polynomial x^5 - 5*x + 2
             sage: s(k.gen())
             a
@@ -1371,10 +1371,10 @@ cdef class LiftingMap(Section):
             sage: F = K.factor(7)[0][0].residue_field()
             sage: L = F.lift_map(); L
             Lifting map:
-              From: Residue field in abar of Fractional ideal (-2*a^4 + a^3 - 4*a^2 + 2*a - 1)
+              From: Residue field in abar of Fractional ideal (2*a^4 - a^3 + 4*a^2 - 2*a + 1)
               To:   Maximal Order in Number Field in a with defining polynomial x^5 + 2
             sage: L.domain()
-            Residue field in abar of Fractional ideal (-2*a^4 + a^3 - 4*a^2 + 2*a - 1)
+            Residue field in abar of Fractional ideal (2*a^4 - a^3 + 4*a^2 - 2*a + 1)
 
             sage: K.<a> = CyclotomicField(7)
             sage: F = K.factor(5)[0][0].residue_field()
@@ -1498,7 +1498,7 @@ cdef class LiftingMap(Section):
             sage: F.<tmod> = K.factor(7)[0][0].residue_field()
             sage: F.lift_map() #indirect doctest
             Lifting map:
-              From: Residue field in tmod of Fractional ideal (-3*theta_12^2 + 1)
+              From: Residue field in tmod of Fractional ideal (theta_12^2 + 2)
               To:   Maximal Order in Cyclotomic Field of order 12 and degree 4
         """
         return "Lifting"
@@ -1515,7 +1515,7 @@ class ResidueFiniteField_prime_modn(ResidueField_generic, FiniteField_prime_modn
         sage: P = K.ideal(29).factor()[1][0]
         sage: k = ResidueField(P)
         sage: k
-        Residue field of Fractional ideal (a^2 + 2*a + 2)
+        Residue field of Fractional ideal (-a^2 - 2*a - 2)
         sage: k.order()
         29
         sage: OK = K.maximal_order()
@@ -1597,7 +1597,7 @@ class ResidueFiniteField_prime_modn(ResidueField_generic, FiniteField_prime_modn
             sage: P = K.ideal(29).factor()[1][0]
             sage: k = ResidueField(P)
             sage: k
-            Residue field of Fractional ideal (a^2 + 2*a + 2)
+            Residue field of Fractional ideal (-a^2 - 2*a - 2)
             sage: OK = K.maximal_order()
             sage: c = OK(a)
             sage: b = k(a); b
diff --git a/src/sage/rings/integer.pyx b/src/sage/rings/integer.pyx
index 77a3a18913..421a90e10e 100644
--- a/src/sage/rings/integer.pyx
+++ b/src/sage/rings/integer.pyx
@@ -5474,7 +5474,7 @@ cdef class Integer(sage.structure.element.EuclideanDomainElement):
             sage: 3._bnfisnorm(QuadraticField(-1, 'i'))
             (1, 3)
             sage: 7._bnfisnorm(CyclotomicField(7))
-            (-zeta7^5 - zeta7^4 - 2*zeta7^3 - zeta7^2 - zeta7 - 1, 1)
+            (zeta7^5 - zeta7^2, 1)
         """
         from sage.rings.rational_field import QQ
         return QQ(self)._bnfisnorm(K, proof=proof, extra_primes=extra_primes)
diff --git a/src/sage/rings/number_field/S_unit_solver.py b/src/sage/rings/number_field/S_unit_solver.py
index 3d6601d536..ba8b5164f0 100644
--- a/src/sage/rings/number_field/S_unit_solver.py
+++ b/src/sage/rings/number_field/S_unit_solver.py
@@ -24,10 +24,10 @@ EXAMPLES::
     sage: from sage.rings.number_field.S_unit_solver import solve_S_unit_equation, eq_up_to_order
     sage: K.<xi> = NumberField(x^2+x+1)
     sage: S = K.primes_above(3)
-    sage: expected = [((2, 1), (4, 0), xi + 2, -xi - 1),
-    ....:             ((5, -1), (4, -1), 1/3*xi + 2/3, -1/3*xi + 1/3),
-    ....:             ((5, 0), (1, 0), -xi, xi + 1),
-    ....:             ((1, 1), (2, 0), -xi + 1, xi)]
+    sage: expected = [((0, 1), (4, 0), xi + 2, -xi - 1),
+    ....:             ((1, -1), (0, -1), 1/3*xi + 2/3, -1/3*xi + 1/3),
+    ....:             ((1, 0), (5, 0), xi + 1, -xi),
+    ....:             ((2, 0), (5, 1), xi, -xi + 1)]
     sage: sols = solve_S_unit_equation(K, S, 200)
     sage: eq_up_to_order(sols, expected)
     True
@@ -1780,20 +1780,20 @@ def sieve_ordering(SUK, q):
         sage: SUK = K.S_unit_group(S=3)
         sage: sieve_data = list(sieve_ordering(SUK, 19))
         sage: sieve_data[0]
-        (Fractional ideal (-2*xi^2 + 3),
-        Fractional ideal (xi - 3),
-        Fractional ideal (2*xi + 1))
+        (Fractional ideal (xi - 3),
+         Fractional ideal (-2*xi^2 + 3),
+         Fractional ideal (2*xi + 1))
 
         sage: sieve_data[1]
-        (Residue field of Fractional ideal (-2*xi^2 + 3),
-        Residue field of Fractional ideal (xi - 3),
-        Residue field of Fractional ideal (2*xi + 1))
+        (Residue field of Fractional ideal (xi - 3),
+         Residue field of Fractional ideal (-2*xi^2 + 3),
+         Residue field of Fractional ideal (2*xi + 1))
 
         sage: sieve_data[2]
-        ([18, 9, 16, 8], [18, 7, 10, 4], [18, 3, 12, 10])
+        ([18, 7, 16, 4], [18, 9, 12, 8], [18, 3, 10, 10])
 
         sage: sieve_data[3]
-        (972, 972, 3888)
+        (486, 648, 11664)
     """
 
     K = SUK.number_field()
@@ -2654,10 +2654,10 @@ def sieve_below_bound(K, S, bound=10, bump=10, split_primes_list=[], verbose=Fal
         sage: S = SUK.primes()
         sage: sols = sieve_below_bound(K, S, 10)
         sage: expected = [
-        ....: ((5, -1), (4, -1), 1/3*xi + 2/3, -1/3*xi + 1/3),
-        ....: ((2, 1), (4, 0), xi + 2, -xi - 1),
-        ....: ((2, 0), (1, 1), xi, -xi + 1),
-        ....: ((5, 0), (1, 0), -xi, xi + 1)]
+        ....: ((1, -1), (0, -1), 1/3*xi + 2/3, -1/3*xi + 1/3),
+        ....: ((0, 1), (4, 0), xi + 2, -xi - 1),
+        ....: ((2, 0), (5, 1), xi, -xi + 1),
+        ....: ((1, 0), (5, 0), xi + 1, -xi)]
         sage: eq_up_to_order(sols, expected)
         True
     """
@@ -2715,10 +2715,10 @@ def solve_S_unit_equation(K, S, prec=106, include_exponents=True, include_bound=
         sage: S = K.primes_above(3)
         sage: sols = solve_S_unit_equation(K, S, 200)
         sage: expected = [
-        ....: ((2, 1), (4, 0), xi + 2, -xi - 1),
-        ....: ((5, -1), (4, -1), 1/3*xi + 2/3, -1/3*xi + 1/3),
-        ....: ((5, 0), (1, 0), -xi, xi + 1),
-        ....: ((1, 1), (2, 0), -xi + 1, xi)]
+        ....: ((0, 1), (4, 0), xi + 2, -xi - 1),
+        ....: ((1, -1), (0, -1), 1/3*xi + 2/3, -1/3*xi + 1/3),
+        ....: ((1, 0), (5, 0), xi + 1, -xi),
+        ....: ((2, 0), (5, 1), xi, -xi + 1)]
         sage: eq_up_to_order(sols, expected)
         True
 
@@ -2726,7 +2726,7 @@ def solve_S_unit_equation(K, S, prec=106, include_exponents=True, include_bound=
 
         sage: solutions, bound = solve_S_unit_equation(K, S, 100, include_bound=True)
         sage: bound
-        6
+        7
 
     You can omit the exponent vectors::
 
diff --git a/src/sage/rings/number_field/class_group.py b/src/sage/rings/number_field/class_group.py
index 46d0ca8c9d..1ad6d583a8 100644
--- a/src/sage/rings/number_field/class_group.py
+++ b/src/sage/rings/number_field/class_group.py
@@ -157,7 +157,7 @@ class FractionalIdealClass(AbelianGroupWithValuesElement):
             sage: C=K.class_group()
             sage: c = C(2, a)
             sage: c^2
-            Fractional ideal class (2, a^2 + 2*a - 1)
+            Fractional ideal class (4, a)
             sage: c^3
             Trivial principal fractional ideal class
             sage: c^1000
@@ -467,7 +467,7 @@ class ClassGroup(AbelianGroupWithValues_class):
             sage: CK = K.class_group()
             sage: CL = L.class_group()
             sage: [CL(I).exponents() for I in CK]
-            [(0,), (4,), (2,)]
+            [(0,), (2,), (4,)]
         """
         if isinstance(args[0], FractionalIdealClass):
             return self.element_class(self, None, self._number_field.ideal(args[0].ideal()))
diff --git a/src/sage/rings/number_field/number_field.py b/src/sage/rings/number_field/number_field.py
index b0fbfffb96..937f7f533d 100644
--- a/src/sage/rings/number_field/number_field.py
+++ b/src/sage/rings/number_field/number_field.py
@@ -3513,7 +3513,7 @@ class NumberField_generic(WithEqualityById, number_field_base.NumberField):
             sage: L.<b> = K.extension(x^2 - 3, x^2 + 1)
             sage: M.<c> = L.extension(x^2 + 1)
             sage: L.ideal(K.ideal(2, a))
-            Fractional ideal (a)
+            Fractional ideal (-a)
             sage: M.ideal(K.ideal(2, a)) == M.ideal(a*(b - c)/2)
             True
 
@@ -4681,7 +4681,7 @@ class NumberField_generic(WithEqualityById, number_field_base.NumberField):
              1/13*a^2 + 7/13*a - 332/13,
              -1/13*a^2 + 6/13*a + 345/13,
              -1,
-             2/13*a^2 + 1/13*a - 755/13]
+             -2/13*a^2 - 1/13*a + 755/13]
             sage: units[5] in (1/13*a^2 - 19/13*a - 7/13, 1/13*a^2 + 20/13*a - 7/13)
             True
             sage: len(units) == 6
@@ -4723,8 +4723,8 @@ class NumberField_generic(WithEqualityById, number_field_base.NumberField):
             sage: K.<a> = QuadraticField(-105)
             sage: K._S_class_group_quotient_matrix((K.ideal(11, a + 4),))
             [0 0]
-            [1 0]
             [0 1]
+            [1 0]
         """
         from sage.matrix.constructor import matrix
         S_clgp_gens = self._S_class_group_and_units(S)[1]
@@ -4845,7 +4845,7 @@ class NumberField_generic(WithEqualityById, number_field_base.NumberField):
              1/13*a^2 + 7/13*a - 332/13,
              -1/13*a^2 + 6/13*a + 345/13,
              -1,
-             2/13*a^2 + 1/13*a - 755/13]
+             -2/13*a^2 - 1/13*a + 755/13]
             sage: gens[5] in (1/13*a^2 - 19/13*a - 7/13, 1/13*a^2 + 20/13*a - 7/13)
             True
             sage: gens[6] in (-1/13*a^2 + 45/13*a - 97/13, 1/13*a^2 - 45/13*a + 97/13)
@@ -6022,28 +6022,37 @@ class NumberField_generic(WithEqualityById, number_field_base.NumberField):
         try:
             return self._integral_basis_dict[v]
         except (AttributeError, KeyError):
-            f = self.pari_polynomial("y")
-            if v:
-                B = f.nfbasis(fa=v)
-            elif self._assume_disc_small:
-                B = f.nfbasis(1)
-            elif not important:
-                # Trial divide the discriminant with primes up to 10^6
-                m = self.pari_polynomial().poldisc().abs().factor(limit=10**6)
-                # Since we only need a *squarefree* factorization for
-                # primes with exponent 1, we need trial division up to D^(1/3)
-                # instead of D^(1/2).
-                trialdivlimit2 = pari(10**12)
-                trialdivlimit3 = pari(10**18)
-                if all(p < trialdivlimit2 or (e == 1 and p < trialdivlimit3) or p.isprime() for p, e in zip(m[0], m[1])):
-                    B = f.nfbasis(fa = m)
-                else:
-                    raise RuntimeError("Unable to factor discriminant with trial division")
+            pass
+
+        f = self.pari_polynomial("y")
+        if v:
+            # NOTE: here we make pari know about potentially big primes factors of
+            # the discriminant, see
+            # https://pari.math.u-bordeaux.fr/cgi-bin/bugreport.cgi?bug=2257
+            primelimit = pari.default("primelimit")
+            primes = [p for p in v if p > primelimit]
+            if primes:
+                pari.addprimes(primes)
+            B = f.nfbasis(fa=v)
+        elif self._assume_disc_small:
+            B = f.nfbasis(1)
+        elif not important:
+            # Trial divide the discriminant with primes up to 10^6
+            m = self.pari_polynomial().poldisc().abs().factor(limit=10**6)
+            # Since we only need a *squarefree* factorization for
+            # primes with exponent 1, we need trial division up to D^(1/3)
+            # instead of D^(1/2).
+            trialdivlimit2 = pari(10**12)
+            trialdivlimit3 = pari(10**18)
+            if all(p < trialdivlimit2 or (e == 1 and p < trialdivlimit3) or p.isprime() for p, e in zip(m[0], m[1])):
+                B = f.nfbasis(fa = m)
             else:
-                B = f.nfbasis()
+                raise RuntimeError("Unable to factor discriminant with trial division")
+        else:
+            B = f.nfbasis()
 
-            self._integral_basis_dict[v] = B
-            return B
+        self._integral_basis_dict[v] = B
+        return B
 
     def reduced_basis(self, prec=None):
         r"""
@@ -6729,7 +6738,8 @@ class NumberField_generic(WithEqualityById, number_field_base.NumberField):
         """
         Return generators for the unit group modulo torsion.
 
-        ALGORITHM: Uses PARI's :pari:`bnfunit` command.
+        ALGORITHM: Uses PARI's :pari:`bnfinit` command (that computes fundamental units
+        among other things).
 
         INPUT:
 
@@ -6751,7 +6761,7 @@ class NumberField_generic(WithEqualityById, number_field_base.NumberField):
             sage: A = x^4 - 10*x^3 + 20*5*x^2 - 15*5^2*x + 11*5^3
             sage: K = NumberField(A, 'a')
             sage: K.units()
-            (8/275*a^3 - 12/55*a^2 + 15/11*a - 3,)
+            (-1/275*a^3 - 4/55*a^2 + 5/11*a - 3,)
 
         For big number fields, provably computing the unit group can
         take a very long time.  In this case, one can ask for the
@@ -6762,14 +6772,14 @@ class NumberField_generic(WithEqualityById, number_field_base.NumberField):
             sage: K.units(proof=True)  # takes forever, not tested
             ...
             sage: K.units(proof=False)  # result not independently verified
-            (a^9 + a - 1,
-             a^15 - a^12 + a^10 - a^9 - 2*a^8 + 3*a^7 + a^6 - 3*a^5 + a^4 + 4*a^3 - 3*a^2 - 2*a + 2,
-             a^16 - a^15 + a^14 - a^12 + a^11 - a^10 - a^8 + a^7 - 2*a^6 + a^4 - 3*a^3 + 2*a^2 - 2*a + 1,
+            (-a^9 - a + 1,
+             -a^16 + a^15 - a^14 + a^12 - a^11 + a^10 + a^8 - a^7 + 2*a^6 - a^4 + 3*a^3 - 2*a^2 + 2*a - 1,
              2*a^16 - a^14 - a^13 + 3*a^12 - 2*a^10 + a^9 + 3*a^8 - 3*a^6 + 3*a^5 + 3*a^4 - 2*a^3 - 2*a^2 + 3*a + 4,
-             2*a^16 - 3*a^15 + 3*a^14 - 3*a^13 + 3*a^12 - a^11 + a^9 - 3*a^8 + 4*a^7 - 5*a^6 + 6*a^5 - 4*a^4 + 3*a^3 - 2*a^2 - 2*a + 4,
-             a^16 - a^15 - 3*a^14 - 4*a^13 - 4*a^12 - 3*a^11 - a^10 + 2*a^9 + 4*a^8 + 5*a^7 + 4*a^6 + 2*a^5 - 2*a^4 - 6*a^3 - 9*a^2 - 9*a - 7,
              a^15 + a^14 + 2*a^11 + a^10 - a^9 + a^8 + 2*a^7 - a^5 + 2*a^3 - a^2 - 3*a + 1,
-             5*a^16 - 6*a^14 + a^13 + 7*a^12 - 2*a^11 - 7*a^10 + 4*a^9 + 7*a^8 - 6*a^7 - 7*a^6 + 8*a^5 + 6*a^4 - 11*a^3 - 5*a^2 + 13*a + 4)
+             -a^16 - a^15 - a^14 - a^13 - a^12 - a^11 - a^10 - a^9 - a^8 - a^7 - a^6 - a^5 - a^4 - a^3 - a^2 + 2,
+             -2*a^16 + 3*a^15 - 3*a^14 + 3*a^13 - 3*a^12 + a^11 - a^9 + 3*a^8 - 4*a^7 + 5*a^6 - 6*a^5 + 4*a^4 - 3*a^3 + 2*a^2 + 2*a - 4,
+             a^15 - a^12 + a^10 - a^9 - 2*a^8 + 3*a^7 + a^6 - 3*a^5 + a^4 + 4*a^3 - 3*a^2 - 2*a + 2,
+             -a^14 - a^13 + a^12 + 2*a^10 + a^8 - 2*a^7 - 2*a^6 + 2*a^3 - a^2 + 2*a - 2)
 
         TESTS:
 
@@ -6778,7 +6788,7 @@ class NumberField_generic(WithEqualityById, number_field_base.NumberField):
 
             sage: K.<a> = NumberField(1/2*x^2 - 1/6)
             sage: K.units()
-            (3*a - 2,)
+            (-3*a + 2,)
         """
         proof = proof_flag(proof)
 
@@ -6796,7 +6806,7 @@ class NumberField_generic(WithEqualityById, number_field_base.NumberField):
                 pass
 
         # get PARI to compute the units
-        B = self.pari_bnf(proof).bnfunit()
+        B = self.pari_bnf(proof).bnf_get_fu()
         B = tuple(self(b, check=False) for b in B)
         if proof:
             # cache the provable results and return them
@@ -6811,7 +6821,7 @@ class NumberField_generic(WithEqualityById, number_field_base.NumberField):
         """
         Return the unit group (including torsion) of this number field.
 
-        ALGORITHM: Uses PARI's :pari:`bnfunit` command.
+        ALGORITHM: Uses PARI's :pari:`bnfinit` and :pari:`bnfunits`.
 
         INPUT:
 
@@ -6857,7 +6867,7 @@ class NumberField_generic(WithEqualityById, number_field_base.NumberField):
             sage: U.gens()
             (u0, u1, u2, u3, u4, u5, u6, u7, u8)
             sage: U.gens_values()  # result not independently verified
-            [-1, a^9 + a - 1, a^15 - a^12 + a^10 - a^9 - 2*a^8 + 3*a^7 + a^6 - 3*a^5 + a^4 + 4*a^3 - 3*a^2 - 2*a + 2, a^16 - a^15 + a^14 - a^12 + a^11 - a^10 - a^8 + a^7 - 2*a^6 + a^4 - 3*a^3 + 2*a^2 - 2*a + 1, 2*a^16 - a^14 - a^13 + 3*a^12 - 2*a^10 + a^9 + 3*a^8 - 3*a^6 + 3*a^5 + 3*a^4 - 2*a^3 - 2*a^2 + 3*a + 4, 2*a^16 - 3*a^15 + 3*a^14 - 3*a^13 + 3*a^12 - a^11 + a^9 - 3*a^8 + 4*a^7 - 5*a^6 + 6*a^5 - 4*a^4 + 3*a^3 - 2*a^2 - 2*a + 4, a^16 - a^15 - 3*a^14 - 4*a^13 - 4*a^12 - 3*a^11 - a^10 + 2*a^9 + 4*a^8 + 5*a^7 + 4*a^6 + 2*a^5 - 2*a^4 - 6*a^3 - 9*a^2 - 9*a - 7, a^15 + a^14 + 2*a^11 + a^10 - a^9 + a^8 + 2*a^7 - a^5 + 2*a^3 - a^2 - 3*a + 1, 5*a^16 - 6*a^14 + a^13 + 7*a^12 - 2*a^11 - 7*a^10 + 4*a^9 + 7*a^8 - 6*a^7 - 7*a^6 + 8*a^5 + 6*a^4 - 11*a^3 - 5*a^2 + 13*a + 4]
+            [-1, -a^9 - a + 1, -a^16 + a^15 - a^14 + a^12 - a^11 + a^10 + a^8 - a^7 + 2*a^6 - a^4 + 3*a^3 - 2*a^2 + 2*a - 1, 2*a^16 - a^14 - a^13 + 3*a^12 - 2*a^10 + a^9 + 3*a^8 - 3*a^6 + 3*a^5 + 3*a^4 - 2*a^3 - 2*a^2 + 3*a + 4, a^15 + a^14 + 2*a^11 + a^10 - a^9 + a^8 + 2*a^7 - a^5 + 2*a^3 - a^2 - 3*a + 1, -a^16 - a^15 - a^14 - a^13 - a^12 - a^11 - a^10 - a^9 - a^8 - a^7 - a^6 - a^5 - a^4 - a^3 - a^2 + 2, -2*a^16 + 3*a^15 - 3*a^14 + 3*a^13 - 3*a^12 + a^11 - a^9 + 3*a^8 - 4*a^7 + 5*a^6 - 6*a^5 + 4*a^4 - 3*a^3 + 2*a^2 + 2*a - 4, a^15 - a^12 + a^10 - a^9 - 2*a^8 + 3*a^7 + a^6 - 3*a^5 + a^4 + 4*a^3 - 3*a^2 - 2*a + 2, -a^14 - a^13 + a^12 + 2*a^10 + a^8 - 2*a^7 - 2*a^6 + 2*a^3 - a^2 + 2*a - 2]
         """
         proof = proof_flag(proof)
 
@@ -7045,7 +7055,7 @@ class NumberField_generic(WithEqualityById, number_field_base.NumberField):
 
             sage: solutions, bound = K.S_unit_solutions(S, prec=100, include_bound=True)
             sage: bound
-            6
+            7
         """
         from .S_unit_solver import solve_S_unit_equation
         return solve_S_unit_equation(self, S, prec, include_exponents, include_bound, proof)
@@ -10884,9 +10894,9 @@ class NumberField_cyclotomic(NumberField_absolute):
         EXAMPLES::
 
             sage: k5.<z> = CyclotomicField(5)
-            sage: gap('E(5)^7 + 3')
+            sage: w = libgap.eval('E(5)^7 + 3')
+            sage: w
             -3*E(5)-2*E(5)^2-3*E(5)^3-3*E(5)^4
-            sage: w = gap('E(5)^7 + 3')
             sage: z^7 + 3
             z^2 + 3
             sage: k5(w) # indirect doctest
@@ -10897,7 +10907,7 @@ class NumberField_cyclotomic(NumberField_absolute):
 
             sage: F = CyclotomicField(8)
             sage: z = F.gen()
-            sage: a = gap(z+1/z); a
+            sage: a = libgap(z+1/z); a
             E(8)-E(8)^3
             sage: F(a)
             -zeta8^3 + zeta8
@@ -10911,6 +10921,7 @@ class NumberField_cyclotomic(NumberField_absolute):
 
         It also works with the old pexpect interface to GAP::
 
+            sage: a = gap(z + 1/z)
             sage: b = gap(Matrix(F,[[z^2,1],[0,a+1]])); b
             [ [ E(4), 1 ], [ 0, 1+E(8)-E(8)^3 ] ]
             sage: b[1,2]
diff --git a/src/sage/rings/number_field/number_field_element.pyx b/src/sage/rings/number_field/number_field_element.pyx
index 0c2499bfee..88bfc6d6d8 100644
--- a/src/sage/rings/number_field/number_field_element.pyx
+++ b/src/sage/rings/number_field/number_field_element.pyx
@@ -1629,7 +1629,7 @@ cdef class NumberFieldElement(FieldElement):
             sage: Q.<X> = K[]
             sage: L.<b> = NumberField(X^4 + a)
             sage: t = (-a).is_norm(L, element=True); t
-            (True, b^3 + 1)
+            (True, -b^3 - 1)
             sage: t[1].norm(K)
             -a
 
@@ -1744,11 +1744,11 @@ cdef class NumberFieldElement(FieldElement):
             sage: Q.<X> = K[]
             sage: L.<b> = NumberField(X^4 + a)
             sage: t = (-a)._rnfisnorm(L); t
-            (b^3 + 1, 1)
+            (-b^3 - 1, 1)
             sage: t[0].norm(K)
             -a
             sage: t = K(3)._rnfisnorm(L); t
-            (-b^3 - a*b^2 - a^2*b + 1, 3*a^2 - 3*a + 6)
+            (b^3 + a*b^2 + a^2*b - 1, 3*a^2 - 3*a + 6)
             sage: t[0].norm(K)*t[1]
             3
 
@@ -1801,7 +1801,7 @@ cdef class NumberFieldElement(FieldElement):
             raise ValueError("L (=%s) must be a relative number field with base field K (=%s) in rnfisnorm" % (L, K))
 
         rnf_data = K.pari_rnfnorm_data(L, proof=proof)
-        x, q = self.__pari__().rnfisnorm(rnf_data)
+        x, q = pari.rnfisnorm(rnf_data, self)
         return L(x, check=False), K(q, check=False)
 
     def _mpfr_(self, R):
diff --git a/src/sage/rings/number_field/number_field_ideal.py b/src/sage/rings/number_field/number_field_ideal.py
index 2e6a368843..aa2cc3c297 100644
--- a/src/sage/rings/number_field/number_field_ideal.py
+++ b/src/sage/rings/number_field/number_field_ideal.py
@@ -231,7 +231,7 @@ class NumberFieldIdeal(Ideal_generic):
             sage: K.<a> = NumberField(x^2 + 3); K
             Number Field in a with defining polynomial x^2 + 3
             sage: f = K.factor(15); f
-            (Fractional ideal (-a))^2 * (Fractional ideal (5))
+            (Fractional ideal (1/2*a + 3/2))^2 * (Fractional ideal (5))
             sage: (f[0][0] < f[1][0])
             True
             sage: (f[0][0] == f[0][0])
@@ -620,7 +620,7 @@ class NumberFieldIdeal(Ideal_generic):
 
             sage: K.<z> = CyclotomicField(7)
             sage: I = K.factor(11)[0][0]; I
-            Fractional ideal (-2*z^4 - 2*z^2 - 2*z + 1)
+            Fractional ideal (-3*z^4 - 2*z^3 - 2*z^2 - 2)
             sage: A = I.free_module()
             sage: A              # warning -- choice of basis can be somewhat random
             Free module of degree 6 and rank 6 over Integer Ring
@@ -3118,7 +3118,7 @@ class NumberFieldFractionalIdeal(MultiplicativeGroupElement, NumberFieldIdeal):
             sage: K.<a> = NumberField(x^5 + 2); K
             Number Field in a with defining polynomial x^5 + 2
             sage: f = K.factor(19); f
-            (Fractional ideal (a^2 + a - 3)) * (Fractional ideal (-2*a^4 - a^2 + 2*a - 1)) * (Fractional ideal (a^2 + a - 1))
+            (Fractional ideal (a^2 + a - 3)) * (Fractional ideal (2*a^4 + a^2 - 2*a + 1)) * (Fractional ideal (a^2 + a - 1))
             sage: [i.residue_class_degree() for i, _ in f]
             [2, 2, 1]
         """
diff --git a/src/sage/rings/number_field/number_field_ideal_rel.py b/src/sage/rings/number_field/number_field_ideal_rel.py
index 37c20150ad..078682fa1f 100644
--- a/src/sage/rings/number_field/number_field_ideal_rel.py
+++ b/src/sage/rings/number_field/number_field_ideal_rel.py
@@ -18,7 +18,7 @@ EXAMPLES::
     sage: G = [from_A(z) for z in I.gens()]; G
     [7, -2*b*a - 1]
     sage: K.fractional_ideal(G)
-    Fractional ideal (2*b*a + 1)
+    Fractional ideal ((1/2*b + 2)*a - 1/2*b + 2)
     sage: K.fractional_ideal(G).absolute_norm().factor()
     7^2
 """
diff --git a/src/sage/rings/number_field/number_field_rel.py b/src/sage/rings/number_field/number_field_rel.py
index 312cbdf874..436af50e87 100644
--- a/src/sage/rings/number_field/number_field_rel.py
+++ b/src/sage/rings/number_field/number_field_rel.py
@@ -396,18 +396,18 @@ class NumberField_relative(NumberField_generic):
             sage: K.<c> = F.extension(Y^2 - (1 + a)*(a + b)*a*b)
             sage: K.subfields(2)
             [
-            (Number Field in c0 with defining polynomial x^2 - 24*x + 72, Ring morphism:
-              From: Number Field in c0 with defining polynomial x^2 - 24*x + 72
+            (Number Field in c0 with defining polynomial x^2 - 24*x + 96, Ring morphism:
+              From: Number Field in c0 with defining polynomial x^2 - 24*x + 96
               To:   Number Field in c with defining polynomial Y^2 + (-2*b - 3)*a - 2*b - 6 over its base field
-              Defn: c0 |--> -6*a + 12, None),
+              Defn: c0 |--> -4*b + 12, None),
             (Number Field in c1 with defining polynomial x^2 - 24*x + 120, Ring morphism:
               From: Number Field in c1 with defining polynomial x^2 - 24*x + 120
               To:   Number Field in c with defining polynomial Y^2 + (-2*b - 3)*a - 2*b - 6 over its base field
               Defn: c1 |--> 2*b*a + 12, None),
-            (Number Field in c2 with defining polynomial x^2 - 24*x + 96, Ring morphism:
-              From: Number Field in c2 with defining polynomial x^2 - 24*x + 96
+            (Number Field in c2 with defining polynomial x^2 - 24*x + 72, Ring morphism:
+              From: Number Field in c2 with defining polynomial x^2 - 24*x + 72
               To:   Number Field in c with defining polynomial Y^2 + (-2*b - 3)*a - 2*b - 6 over its base field
-              Defn: c2 |--> -4*b + 12, None)
+              Defn: c2 |--> -6*a + 12, None)
             ]
             sage: K.subfields(8, 'w')
             [
diff --git a/src/sage/rings/number_field/order.py b/src/sage/rings/number_field/order.py
index b5cdeb2a96..be1548f7a2 100644
--- a/src/sage/rings/number_field/order.py
+++ b/src/sage/rings/number_field/order.py
@@ -2180,7 +2180,7 @@ def EisensteinIntegers(names="omega"):
         sage: R
         Eisenstein Integers in Number Field in omega with defining polynomial x^2 + x + 1 with omega = -0.50000000000000000? + 0.866025403784439?*I
         sage: factor(3 + omega)
-        (omega) * (-3*omega - 2)
+        (-1) * (-omega - 3)
         sage: CC(omega)
         -0.500000000000000 + 0.866025403784439*I
         sage: omega.minpoly()
diff --git a/src/sage/rings/number_field/unit_group.py b/src/sage/rings/number_field/unit_group.py
index 5c94f7407d..2a3c617ca3 100644
--- a/src/sage/rings/number_field/unit_group.py
+++ b/src/sage/rings/number_field/unit_group.py
@@ -15,12 +15,12 @@ The first generator is a primitive root of unity in the field::
     sage: UK.gens_values()  # random
     [-1/12*a^3 + 1/6*a, 1/24*a^3 + 1/4*a^2 - 1/12*a - 1]
     sage: UK.gen(0).value()
-    -1/12*a^3 + 1/6*a
+    1/12*a^3 - 1/6*a
 
     sage: UK.gen(0)
     u0
     sage: UK.gen(0) + K.one()   # coerce abstract generator into number field
-    -1/12*a^3 + 1/6*a + 1
+    1/12*a^3 - 1/6*a + 1
 
     sage: [u.multiplicative_order() for u in UK.gens()]
     [4, +Infinity]
@@ -37,18 +37,18 @@ as elements of an abstract multiplicative group::
     sage: UK(-1)
     u0^2
     sage: [UK(u) for u in (x^4-1).roots(K, multiplicities=False)]
-    [1, u0^2, u0^3, u0]
+    [1, u0^2, u0, u0^3]
 
     sage: UK.fundamental_units() # random
     [1/24*a^3 + 1/4*a^2 - 1/12*a - 1]
     sage: torsion_gen = UK.torsion_generator();  torsion_gen
     u0
     sage: torsion_gen.value()
-    -1/12*a^3 + 1/6*a
+    1/12*a^3 - 1/6*a
     sage: UK.zeta_order()
     4
     sage: UK.roots_of_unity()
-    [-1/12*a^3 + 1/6*a, -1, 1/12*a^3 - 1/6*a, 1]
+    [1/12*a^3 - 1/6*a, -1, -1/12*a^3 + 1/6*a, 1]
 
 Exp and log functions provide maps between units as field elements and exponent
 vectors with respect to the generators::
@@ -82,7 +82,7 @@ S-unit groups may be constructed, where S is a set of primes::
     sage: SUK.rank()
     4
     sage: SUK.gens_values()
-    [-1, a^2 + 1, a^5 + a^4 - a^2 - a - 1, a + 1, -a + 1]
+    [-1, a^2 + 1, -a^5 - a^4 + a^2 + a + 1, a + 1, a - 1]
     sage: u = 9*prod(SUK.gens_values()); u
     -18*a^5 - 18*a^4 - 18*a^3 - 9*a^2 + 9*a + 27
     sage: SUK.log(u)
@@ -100,29 +100,29 @@ A relative number field example::
     sage: UL.zeta_order()
     24
     sage: UL.roots_of_unity()
-    [-b*a - b,
-     b^2*a,
-     b^3,
-     a + 1,
-     -b*a,
-     -b^2,
-     b^3*a + b^3,
-     a,
-     b,
+    [-b*a,
      -b^2*a - b^2,
-     b^3*a,
-     -1,
-     b*a + b,
-     -b^2*a,
      -b^3,
-     -a - 1,
-     b*a,
-     b^2,
-     -b^3*a - b^3,
      -a,
+     -b*a - b,
+     -b^2,
+     b^3*a,
+     -a - 1,
      -b,
+     b^2*a,
+     b^3*a + b^3,
+     -1,
+     b*a,
      b^2*a + b^2,
+     b^3,
+     a,
+     b*a + b,
+     b^2,
      -b^3*a,
+     a + 1,
+     b,
+     -b^2*a,
+     -b^3*a - b^3,
      1]
 
 A relative extension example, which worked thanks to the code review by F.W.Clarke::
@@ -199,7 +199,7 @@ class UnitGroup(AbelianGroupWithValues_class):
         sage: UK.gen(5)
         u5
         sage: UK.gen(5).value()
-        z^7 + z
+        -z^7 - z
 
     An S-unit group::
 
@@ -216,7 +216,7 @@ class UnitGroup(AbelianGroupWithValues_class):
         sage: SUK.zeta_order()
         26
         sage: SUK.log(21*z)
-        (12, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1)
+        (25, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1)
     """
     # This structure is not a parent in the usual sense. The
     # "elements" are NumberFieldElement_absolute. Instead, they should
@@ -250,7 +250,7 @@ class UnitGroup(AbelianGroupWithValues_class):
             sage: UK.gens()
             (u0, u1)
             sage: UK.gens_values()
-            [-1, 6*a - 37]
+            [-1, -6*a + 37]
 
             sage: K.<a> = QuadraticField(-3)
             sage: UK = K.unit_group(); UK
@@ -258,7 +258,7 @@ class UnitGroup(AbelianGroupWithValues_class):
             sage: UK.gens()
             (u,)
             sage: UK.gens_values()
-            [1/2*a + 1/2]
+            [-1/2*a + 1/2]
 
             sage: K.<z> = CyclotomicField(13)
             sage: UK = K.unit_group(); UK
@@ -323,34 +323,31 @@ class UnitGroup(AbelianGroupWithValues_class):
             self.__S = S
             self.__pS = pS = [P.pari_prime() for P in S]
 
-        # compute the fundamental units via pari:
-        fu = [K(u, check=False) for u in pK.bnfunit()]
-        self.__nfu = len(fu)
-
-        # compute the additional S-unit generators:
+        # compute units
+        # NOTE: old pari syntax for S-units (< 2.13.0): pK.bnfsunit(pS)
+        # NOTE: pari >= 2.13.0: the first component of the result of bnfunits
+        # are *all* units starting with S-units, followed by fundamental units
+        # followed by the torsion unit.
         if S:
-            self.__S_unit_data = pK.bnfsunit(pS)
-            su = [K(u, check=False) for u in self.__S_unit_data[0]]
+            self.__S_unit_data = pK.bnfunits(pS)
         else:
-            su = []
-        self.__nsu = len(su)
+            self.__S_unit_data = pK.bnfunits()
+        # TODO: converting the factored matrix representation of bnfunits into polynomial
+        # form is a *big* waste of time
+        su_fu_tu = [pK.nfbasistoalg(pK.nffactorback(z)) for z in self.__S_unit_data[0]]
+
+        self.__nfu = len(pK.bnf_get_fu())           # number of fundamental units
+        self.__nsu = len(su_fu_tu) - self.__nfu - 1 # number of S-units
+        self.__ntu = pK.bnf_get_tu()[0]             # order of torsion
         self.__rank = self.__nfu + self.__nsu
 
-        # compute a torsion generator and pick the 'simplest' one:
-        n, z = pK[7][3] # number of roots of unity and bnf.tu as in pari documentation
-        n = ZZ(n)
-        self.__ntu = n
-        z = K(z, check=False)
-
-        # If we replaced z by another torsion generator we would need
-        # to allow for this in the dlog function!  So we do not.
+        # Move the torsion unit first, then fundamental units then S-units
+        gens = [K(u, check=False) for u in su_fu_tu]
+        gens = [gens[-1]] + gens[self.__nsu:-1] + gens[:self.__nsu]
 
-        # Store the actual generators (torsion first):
-        gens = [z] + fu + su
-        values = Sequence(gens, immutable=True, universe=self, check=False)
         # Construct the abtract group:
-        gens_orders = tuple([ZZ(n)]+[ZZ(0)]*(self.__rank))
-        AbelianGroupWithValues_class.__init__(self, gens_orders, 'u', values, number_field)
+        gens_orders = tuple([ZZ(self.__ntu)]+[ZZ(0)]*(self.__rank))
+        AbelianGroupWithValues_class.__init__(self, gens_orders, 'u', gens, number_field)
 
     def _element_constructor_(self, u):
         """
@@ -375,7 +372,7 @@ class UnitGroup(AbelianGroupWithValues_class):
             sage: UK.gens()
             (u0, u1)
             sage: UK.gens_values()
-            [-1, 6*a - 37]
+            [-1, -6*a + 37]
             sage: UK.ngens()
             2
             sage: [UK(u) for u in UK.gens()]
@@ -394,8 +391,8 @@ class UnitGroup(AbelianGroupWithValues_class):
         except TypeError:
             raise ValueError("%s is not an element of %s"%(u,K))
         if self.__S:
-            m = pK.bnfissunit(self.__S_unit_data, pari(u)).mattranspose()
-            if m.ncols()==0:
+            m = pK.bnfisunit(pari(u), self.__S_unit_data).mattranspose()
+            if m.ncols() == 0:
                 raise ValueError("%s is not an S-unit"%u)
         else:
             if not u.is_integral() or u.norm().abs() != 1:
@@ -405,9 +402,8 @@ class UnitGroup(AbelianGroupWithValues_class):
         # convert column matrix to a list:
         m = [ZZ(m[0,i].sage()) for i in range(m.ncols())]
 
-        # NB pari puts the torsion after the fundamental units, before
-        # the extra S-units but we have the torsion first:
-        m = [m[self.__nfu]] + m[:self.__nfu] + m[self.__nfu+1:]
+        # NOTE: pari ordering for the units is (S-units, fundamental units, torsion unit)
+        m = [m[-1]] + m[self.__nsu:-1] + m[:self.__nsu]
 
         return self.element_class(self, m)
 
@@ -527,9 +523,9 @@ class UnitGroup(AbelianGroupWithValues_class):
             sage: U.zeta(2, all=True)
             [-1]
             sage: U.zeta(3)
-            1/2*z - 1/2
+            -1/2*z - 1/2
             sage: U.zeta(3, all=True)
-            [1/2*z - 1/2, -1/2*z - 1/2]
+            [-1/2*z - 1/2, 1/2*z - 1/2]
             sage: U.zeta(4)
             Traceback (most recent call last):
             ...
@@ -645,7 +641,7 @@ class UnitGroup(AbelianGroupWithValues_class):
            sage: SUK = UnitGroup(K,S=2)
            sage: v = (3,1,4,1,5,9,2)
            sage: u = SUK.exp(v); u
-           -8732*z^11 + 15496*z^10 + 51840*z^9 + 68804*z^8 + 51840*z^7 + 15496*z^6 - 8732*z^5 + 34216*z^3 + 64312*z^2 + 64312*z + 34216
+           8732*z^11 - 15496*z^10 - 51840*z^9 - 68804*z^8 - 51840*z^7 - 15496*z^6 + 8732*z^5 - 34216*z^3 - 64312*z^2 - 64312*z - 34216
            sage: SUK.log(u)
            (3, 1, 4, 1, 5, 9, 2)
            sage: SUK.log(u) == v
@@ -692,7 +688,7 @@ class UnitGroup(AbelianGroupWithValues_class):
            sage: SUK = UnitGroup(K,S=2)
            sage: v = (3,1,4,1,5,9,2)
            sage: u = SUK.exp(v); u
-           -8732*z^11 + 15496*z^10 + 51840*z^9 + 68804*z^8 + 51840*z^7 + 15496*z^6 - 8732*z^5 + 34216*z^3 + 64312*z^2 + 64312*z + 34216
+           8732*z^11 - 15496*z^10 - 51840*z^9 - 68804*z^8 - 51840*z^7 - 15496*z^6 + 8732*z^5 - 34216*z^3 - 64312*z^2 - 64312*z - 34216
            sage: SUK.log(u)
            (3, 1, 4, 1, 5, 9, 2)
            sage: SUK.log(u) == v
diff --git a/src/sage/rings/polynomial/multi_polynomial_ring.py b/src/sage/rings/polynomial/multi_polynomial_ring.py
index 077f380730..ce6e517c78 100644
--- a/src/sage/rings/polynomial/multi_polynomial_ring.py
+++ b/src/sage/rings/polynomial/multi_polynomial_ring.py
@@ -525,6 +525,8 @@ class MPolynomialRing_polydict( MPolynomialRing_macaulay2_repr, PolynomialRing_s
             # univariate polynomials.  Below, v is the variable
             # with highest priority, and the x[i] are expressions
             # in the remaining variables.
+            if x == 0:
+                return self.zero()
             v = self.gens_dict_recursive()[str(x.variable())]
             return sum(self(x[i]) * v ** i for i in range(x.poldegree() + 1))
 
diff --git a/src/sage/rings/polynomial/polynomial_element.pyx b/src/sage/rings/polynomial/polynomial_element.pyx
index b7f3058ce7..4615107773 100644
--- a/src/sage/rings/polynomial/polynomial_element.pyx
+++ b/src/sage/rings/polynomial/polynomial_element.pyx
@@ -7550,7 +7550,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
             [(-3.5074662110434039?e451, 1)]
             sage: p = bigc*x + 1
             sage: p.roots(ring=RR)
-            [(0.000000000000000, 1)]
+            [(-2.85106096489671e-452, 1)]
             sage: p.roots(ring=AA)
             [(-2.8510609648967059?e-452, 1)]
             sage: p.roots(ring=QQbar)
diff --git a/src/sage/rings/polynomial/polynomial_quotient_ring.py b/src/sage/rings/polynomial/polynomial_quotient_ring.py
index 1ed1faaf15..6d4bb74464 100644
--- a/src/sage/rings/polynomial/polynomial_quotient_ring.py
+++ b/src/sage/rings/polynomial/polynomial_quotient_ring.py
@@ -1305,7 +1305,16 @@ class PolynomialQuotientRing_generic(CommutativeRing):
         fixed in :trac:`14489`)::
 
             sage: S.S_class_group([K.ideal(a)])
-            [((1/4*xbar^2 + 31/4, (-1/8*a + 1/8)*xbar^2 - 31/8*a + 31/8, 1/16*xbar^3 + 1/16*xbar^2 + 31/16*xbar + 31/16, -1/16*a*xbar^3 + (1/16*a + 1/8)*xbar^2 - 31/16*a*xbar + 31/16*a + 31/8), 6), ((-1/4*xbar^2 - 23/4, (1/8*a - 1/8)*xbar^2 + 23/8*a - 23/8, -1/16*xbar^3 - 1/16*xbar^2 - 23/16*xbar - 23/16, 1/16*a*xbar^3 + (-1/16*a - 1/8)*xbar^2 + 23/16*a*xbar - 23/16*a - 23/8), 2)]
+            [((1/4*xbar^2 + 31/4,
+               (-1/8*a + 1/8)*xbar^2 - 31/8*a + 31/8,
+               1/16*xbar^3 + 1/16*xbar^2 + 31/16*xbar + 31/16,
+               -1/16*a*xbar^3 + (1/16*a + 1/8)*xbar^2 - 31/16*a*xbar + 31/16*a + 31/8),
+              6),
+             ((-1/4*xbar^2 - 23/4,
+               1/4*a*xbar^2 + 23/4*a,
+               -1/8*xbar^3 + 1/8*xbar^2 - 23/8*xbar + 23/8,
+               (1/16*a - 1/16)*xbar^3 + (-1/16*a + 1/16)*xbar^2 + (23/16*a - 23/16)*xbar - 23/16*a + 23/16),
+              2)]
 
         Note that all the returned values live where we expect them to::
 
diff --git a/src/sage/schemes/elliptic_curves/ell_finite_field.py b/src/sage/schemes/elliptic_curves/ell_finite_field.py
index 87fbec3b69..adf0d2c73f 100644
--- a/src/sage/schemes/elliptic_curves/ell_finite_field.py
+++ b/src/sage/schemes/elliptic_curves/ell_finite_field.py
@@ -785,11 +785,11 @@ class EllipticCurve_finite_field(EllipticCurve_field, HyperellipticCurve_finite_
             sage: len(E.gens())
             2
             sage: E.cardinality()
-            867361737988403547207212930746733987710588
+            867361737988403547206134229616487867594472
             sage: E.gens()[0].order()
-            433680868994201773603606465373366993855294
+            433680868994201773603067114808243933797236
             sage: E.gens()[1].order()
-            433680868994201773603606465373366993855294
+            433680868994201773603067114808243933797236
         """
         G = self.__pari__().ellgroup(flag=1)
         return tuple(self.point(list(pt)) for pt in G[2])
diff --git a/src/sage/schemes/elliptic_curves/ell_generic.py b/src/sage/schemes/elliptic_curves/ell_generic.py
index ac5b4a98b4..194501bc94 100644
--- a/src/sage/schemes/elliptic_curves/ell_generic.py
+++ b/src/sage/schemes/elliptic_curves/ell_generic.py
@@ -2934,15 +2934,7 @@ class EllipticCurve_generic(WithEqualityById, plane_curve.ProjectivePlaneCurve):
             sage: K.<a> = QuadraticField(2)
             sage: E = EllipticCurve([1,a])
             sage: E.pari_curve()
-            [Mod(0, y^2 - 2), Mod(0, y^2 - 2), Mod(0, y^2 - 2), Mod(1, y^2 - 2),
-            Mod(y, y^2 - 2), Mod(0, y^2 - 2), Mod(2, y^2 - 2), Mod(4*y, y^2 - 2),
-            Mod(-1, y^2 - 2), Mod(-48, y^2 - 2), Mod(-864*y, y^2 - 2),
-            Mod(-928, y^2 - 2), Mod(3456/29, y^2 - 2), Vecsmall([5]),
-            [[y^2 - 2, [2, 0], 8, 1, [[1, -1.41421356237310;
-            1, 1.41421356237310], [1, -1.41421356237310; 1, 1.41421356237310],
-            [1, -1; 1, 1], [2, 0; 0, 4], [4, 0; 0, 2], [2, 0; 0, 1],
-            [2, [0, 2; 1, 0]], []], [-1.41421356237310, 1.41421356237310],
-            [1, y], [1, 0; 0, 1], [1, 0, 0, 2; 0, 1, 1, 0]]], [0, 0, 0, 0, 0]]
+            [Mod(0, y^2 - 2), Mod(0, y^2 - 2), Mod(0, y^2 - 2), Mod(1, y^2 - 2), Mod(y, y^2 - 2), Mod(0, y^2 - 2), Mod(2, y^2 - 2), Mod(4*y, y^2 - 2), Mod(-1, y^2 - 2), Mod(-48, y^2 - 2), Mod(-864*y, y^2 - 2), Mod(-928, y^2 - 2), Mod(3456/29, y^2 - 2), Vecsmall([5]), [[y^2 - 2, [2, 0], 8, 1, [[1, -1.41421356237310; 1, 1.41421356237310], [1, -1.41421356237310; 1, 1.41421356237310], [16, -23; 16, 23], [2, 0; 0, 4], [4, 0; 0, 2], [2, 0; 0, 1], [2, [0, 2; 1, 0]], [2]], [-1.41421356237310, 1.41421356237310], [1, y], [1, 0; 0, 1], [1, 0, 0, 2; 0, 1, 1, 0]]], [0, 0, 0, 0, 0]]
 
         PARI no longer requires that the `j`-invariant has negative `p`-adic valuation::
 
diff --git a/src/sage/schemes/elliptic_curves/ell_number_field.py b/src/sage/schemes/elliptic_curves/ell_number_field.py
index dded480ba8..26f2d9c9d4 100644
--- a/src/sage/schemes/elliptic_curves/ell_number_field.py
+++ b/src/sage/schemes/elliptic_curves/ell_number_field.py
@@ -214,9 +214,9 @@ class EllipticCurve_number_field(EllipticCurve_field):
             sage: E == loads(dumps(E))
             True
             sage: E.simon_two_descent()
-            (2, 2, [(0 : 0 : 1), (1/8*a + 5/8 : -3/16*a - 7/16 : 1)])
+            (2, 2, [(0 : 0 : 1)])
             sage: E.simon_two_descent(lim1=3, lim3=20, limtriv=5, maxprob=7, limbigprime=10)
-            (2, 2, [(-1 : 0 : 1), (-1/8*a + 5/8 : -3/16*a - 9/16 : 1)])
+            (2, 2, [(-1 : 0 : 1), (5/32*a + 1/32 : -19/128*a + 41/128 : 1)])
 
         ::
 
@@ -244,22 +244,22 @@ class EllipticCurve_number_field(EllipticCurve_field):
               C = Mod(y, y^2 + 7)
             <BLANKLINE>
               Computing L(S,2)
-              L(S,2) = [Mod(Mod(-1/2*y + 1/2, y^2 + 7)*x^2 + Mod(-1/2*y - 1/2, y^2 + 7)*x + Mod(-y - 1, y^2 + 7), x^3 + Mod(1, y^2 + 7)*x + Mod(y, y^2 + 7)), Mod(-x^2 + Mod(-1/2*y - 1/2, y^2 + 7)*x + 1, x^3 + Mod(1, y^2 + 7)*x + Mod(y, y^2 + 7)), Mod(-1, x^3 + Mod(1, y^2 + 7)*x + Mod(y, y^2 + 7)), Mod(x^2 + 2, x^3 + Mod(1, y^2 + 7)*x + Mod(y, y^2 + 7)), Mod(x + Mod(1/2*y + 3/2, y^2 + 7), x^3 + Mod(1, y^2 + 7)*x + Mod(y, y^2 + 7)), Mod(x + Mod(1/2*y - 3/2, y^2 + 7), x^3 + Mod(1, y^2 + 7)*x + Mod(y, y^2 + 7))]
+              L(S,2) = [Mod(Mod(1/2*y - 1/2, y^2 + 7)*x^2 + Mod(1/2*y + 1/2, y^2 + 7)*x + Mod(y + 1, y^2 + 7), x^3 + Mod(1, y^2 + 7)*x + Mod(y, y^2 + 7)), Mod(x^2 + Mod(1/2*y + 1/2, y^2 + 7)*x - 1, x^3 + Mod(1, y^2 + 7)*x + Mod(y, y^2 + 7)), Mod(-1, x^3 + Mod(1, y^2 + 7)*x + Mod(y, y^2 + 7)), Mod(x^2 + 2, x^3 + Mod(1, y^2 + 7)*x + Mod(y, y^2 + 7)), Mod(x + Mod(1/2*y + 3/2, y^2 + 7), x^3 + Mod(1, y^2 + 7)*x + Mod(y, y^2 + 7)), Mod(x + Mod(1/2*y - 3/2, y^2 + 7), x^3 + Mod(1, y^2 + 7)*x + Mod(y, y^2 + 7))]
             <BLANKLINE>
               Computing the Selmer group
               #LS2gen = 2
-               LS2gen = [Mod(Mod(-1/2*y + 1/2, y^2 + 7)*x^2 + Mod(-1/2*y - 1/2, y^2 + 7)*x + Mod(-y - 1, y^2 + 7), x^3 + Mod(1, y^2 + 7)*x + Mod(y, y^2 + 7)), Mod(x^2 + Mod(1/2*y + 1/2, y^2 + 7)*x - 1, x^3 + Mod(1, y^2 + 7)*x + Mod(y, y^2 + 7))]
+               LS2gen = [Mod(x^2 + Mod(1/2*y + 1/2, y^2 + 7)*x - 1, x^3 + Mod(1, y^2 + 7)*x + Mod(y, y^2 + 7)), Mod(Mod(-1/2*y + 1/2, y^2 + 7)*x^2 + Mod(-1/2*y - 1/2, y^2 + 7)*x + Mod(-y - 1, y^2 + 7), x^3 + Mod(1, y^2 + 7)*x + Mod(y, y^2 + 7))]
               Search for trivial points on the curve
              Trivial points on the curve = [[Mod(1/2*y + 3/2, y^2 + 7), Mod(-y - 2, y^2 + 7)], [1, 1, 0], [Mod(1/2*y + 3/2, y^2 + 7), Mod(-y - 2, y^2 + 7), 1]]
+              zc = Mod(Mod(1, y^2 + 7)*x^2 + Mod(1/2*y + 1/2, y^2 + 7)*x + Mod(-1, y^2 + 7), x^3 + Mod(1, y^2 + 7)*x + Mod(y, y^2 + 7))
+              comes from the trivial point [Mod(1/2*y + 3/2, y^2 + 7), Mod(-y - 2, y^2 + 7)]
               zc = Mod(Mod(-1/2*y + 1/2, y^2 + 7)*x^2 + Mod(-1/2*y - 1/2, y^2 + 7)*x + Mod(-y - 1, y^2 + 7), x^3 + Mod(1, y^2 + 7)*x + Mod(y, y^2 + 7))
-              Hilbert symbol (Mod(1, y^2 + 7),Mod(-2*y + 2, y^2 + 7)) =
+              Hilbert symbol (Mod(1, y^2 + 7),Mod(-2*y + 2, y^2 + 7)) = 
               sol of quadratic equation = [1, 1, 0]~
               zc*z1^2 = Mod(4*x + Mod(-2*y + 6, y^2 + 7), x^3 + Mod(1, y^2 + 7)*x + Mod(y, y^2 + 7))
               quartic: (-1)*Y^2 = x^4 + (3*y - 9)*x^2 + (-8*y + 16)*x + (9/2*y - 11/2)
               reduced: Y^2 = -x^4 + (-3*y + 9)*x^2 + (-8*y + 16)*x + (-9/2*y + 11/2)
               not ELS at [2, [0, 1]~, 1, 1, [1, -2; 1, 0]]
-              zc = Mod(Mod(1, y^2 + 7)*x^2 + Mod(1/2*y + 1/2, y^2 + 7)*x + Mod(-1, y^2 + 7), x^3 + Mod(1, y^2 + 7)*x + Mod(y, y^2 + 7))
-              comes from the trivial point [Mod(1/2*y + 3/2, y^2 + 7), Mod(-y - 2, y^2 + 7)]
               m1 = 1
               m2 = 1
             #S(E/K)[2]    = 2
@@ -298,8 +298,8 @@ class EllipticCurve_number_field(EllipticCurve_field):
             sage: E.simon_two_descent()  # long time (4s on sage.math, 2013)
             (3,
              3,
-             [(0 : 0 : 1),
-              (-1/2*zeta43_0^2 - 1/2*zeta43_0 + 7 : -3/2*zeta43_0^2 - 5/2*zeta43_0 + 18 : 1)...)
+             [(5/8*zeta43_0^2 + 17/8*zeta43_0 - 9/4 : -27/16*zeta43_0^2 - 103/16*zeta43_0 + 39/8 : 1),
+              (0 : 0 : 1)])
         """
         verbose = int(verbose)
         if known_points is None:
diff --git a/src/sage/schemes/elliptic_curves/ell_rational_field.py b/src/sage/schemes/elliptic_curves/ell_rational_field.py
index 2953e80cac..b9d34cd64a 100644
--- a/src/sage/schemes/elliptic_curves/ell_rational_field.py
+++ b/src/sage/schemes/elliptic_curves/ell_rational_field.py
@@ -5407,9 +5407,9 @@ class EllipticCurve_rational_field(EllipticCurve_number_field):
 
             sage: E = EllipticCurve('37a1')
             sage: E.eval_modular_form([1.5+I,2.0+I,2.5+I],100) # abs tol 1e-20
-            [-0.0018743978548152085771342944989052703431,
-             0.0018604485340371083710285594393397945456,
-            -0.0018743978548152085771342944989052703431]
+            [-0.00187439785481520858 - 6.91083670607514589e-22*I,
+             0.00186044853403710837 + 3.71914507780688601e-22*I,
+             -0.00187439785481520858 - 6.39417173217386647e-23*I]
 
             sage: E.eval_modular_form(2.1+I, 100) # abs tol 1e-16
             [0.00150864362757267079 + 0.00109100341113449845*I]
diff --git a/src/sage/schemes/elliptic_curves/isogeny_small_degree.py b/src/sage/schemes/elliptic_curves/isogeny_small_degree.py
index 7178da3685..ea0a7e88bc 100644
--- a/src/sage/schemes/elliptic_curves/isogeny_small_degree.py
+++ b/src/sage/schemes/elliptic_curves/isogeny_small_degree.py
@@ -801,8 +801,8 @@ def isogenies_5_0(E, minimal_models=True):
         sage: K.<a> = NumberField(x**6-320*x**3-320)
         sage: E = EllipticCurve(K,[0,0,1,0,0])
         sage: isogenies_5_0(E)
-        [Isogeny of degree 5 from Elliptic Curve defined by y^2 + y = x^3 over Number Field in a with defining polynomial x^6 - 320*x^3 - 320 to Elliptic Curve defined by y^2 + y = x^3 + (643/8*a^5-15779/48*a^4-32939/24*a^3-71989/2*a^2+214321/6*a-112115/3)*x + (2901961/96*a^5+4045805/48*a^4+12594215/18*a^3-30029635/6*a^2+15341626/3*a-38944312/9) over Number Field in a with defining polynomial x^6 - 320*x^3 - 320,
-        Isogeny of degree 5 from Elliptic Curve defined by y^2 + y = x^3 over Number Field in a with defining polynomial x^6 - 320*x^3 - 320 to Elliptic Curve defined by y^2 + y = x^3 + (-1109/8*a^5-53873/48*a^4-180281/24*a^3-14491/2*a^2+35899/6*a-43745/3)*x + (-17790679/96*a^5-60439571/48*a^4-77680504/9*a^3+1286245/6*a^2-4961854/3*a-73854632/9) over Number Field in a with defining polynomial x^6 - 320*x^3 - 320]
+        [Isogeny of degree 5 from Elliptic Curve defined by y^2 + y = x^3 over Number Field in a with defining polynomial x^6 - 320*x^3 - 320 to Elliptic Curve defined by y^2 + y = x^3 + (241565/32*a^5-362149/48*a^4+180281/24*a^3-9693307/4*a^2+14524871/6*a-7254985/3)*x + (1660391123/192*a^5-829315373/96*a^4+77680504/9*a^3-66622345345/24*a^2+33276655441/12*a-24931615912/9) over Number Field in a with defining polynomial x^6 - 320*x^3 - 320,
+         Isogeny of degree 5 from Elliptic Curve defined by y^2 + y = x^3 over Number Field in a with defining polynomial x^6 - 320*x^3 - 320 to Elliptic Curve defined by y^2 + y = x^3 + (47519/32*a^5-72103/48*a^4+32939/24*a^3-1909753/4*a^2+2861549/6*a-1429675/3)*x + (-131678717/192*a^5+65520419/96*a^4-12594215/18*a^3+5280985135/24*a^2-2637787519/12*a+1976130088/9) over Number Field in a with defining polynomial x^6 - 320*x^3 - 320]
     """
     F = E.base_field()
     if E.j_invariant() != 0:
diff --git a/src/sage/schemes/elliptic_curves/period_lattice.py b/src/sage/schemes/elliptic_curves/period_lattice.py
index 42793dc1b3..6ae1dd593d 100644
--- a/src/sage/schemes/elliptic_curves/period_lattice.py
+++ b/src/sage/schemes/elliptic_curves/period_lattice.py
@@ -1626,7 +1626,7 @@ class PeriodLattice_ell(PeriodLattice):
             sage: embs = K.embeddings(CC)
             sage: Lambda = E.period_lattice(embs[0])
             sage: Lambda.elliptic_logarithm(P+3*Q, 100)
-            4.7100131126199672766973600998
+            4.3543876242043418255250464574
             sage: R.<x> = QQ[]
             sage: K.<a> = NumberField(x^2 + x + 5)
             sage: E = EllipticCurve(K, [0,0,1,-3,-5])
diff --git a/src/sage/schemes/toric/chow_group.py b/src/sage/schemes/toric/chow_group.py
index 89e82467b0..78730f930a 100644
--- a/src/sage/schemes/toric/chow_group.py
+++ b/src/sage/schemes/toric/chow_group.py
@@ -62,7 +62,7 @@ EXAMPLES::
     7
     sage: a = sum( A.gen(i) * (i+1) for i in range(A.ngens()) )   # an element of A
     sage: a  # long time (2s on sage.math, 2011)
-    ( 3 | 1 mod 7 | 0 mod 2, 1 mod 2, 4, 5, 6, 7, 8 | 9 )
+    ( 3 | 1 mod 7 | 1 mod 2, 0 mod 2, 4, 5, 6, 7, 8 | 9 )
 
 The Chow group elements are printed as ``( a0 | a1 mod 7 | a2 mod 2,
 a3 mod 2, a4, a5, a6, a7, a8 | a9 )``, which denotes the element of
@@ -93,13 +93,13 @@ Cones of toric varieties can determine their own Chow cycle::
     sage: cone = X.fan(dim=2)[3]; cone
     2-d cone of Rational polyhedral fan in 3-d lattice N
     sage: A_cone = A(cone); A_cone
-    ( 0 | 1 mod 7 | 0 mod 2, 0 mod 2, 0, 0, 0, 0, 0 | 0 )
+    ( 0 | 6 mod 7 | 0 mod 2, 0 mod 2, 0, 0, 0, 0, 0 | 0 )
     sage: A_cone.degree()
     1
     sage: 2 * A_cone
-    ( 0 | 2 mod 7 | 0 mod 2, 0 mod 2, 0, 0, 0, 0, 0 | 0 )
+    ( 0 | 5 mod 7 | 0 mod 2, 0 mod 2, 0, 0, 0, 0, 0 | 0 )
     sage: A_cone + A.gen(0)
-    ( 0 | 1 mod 7 | 0 mod 2, 1 mod 2, 0, 0, 0, 0, 0 | 0 )
+    ( 0 | 6 mod 7 | 1 mod 2, 0 mod 2, 0, 0, 0, 0, 0 | 0 )
 
 Chow cycles can be of mixed degrees::
 
diff --git a/src/sage/schemes/toric/homset.py b/src/sage/schemes/toric/homset.py
index 4bff92bcb0..ff6211cabe 100644
--- a/src/sage/schemes/toric/homset.py
+++ b/src/sage/schemes/toric/homset.py
@@ -467,12 +467,27 @@ class SchemeHomset_points_toric_field(SchemeHomset_points_toric_base):
         sage: point_set.cardinality()
         21
         sage: sorted(X.point_set().list())
-        [[0 : 0 : 1], [0 : 1 : 0], [0 : 1 : 1], [0 : 1 : 3],
-         [1 : 0 : 0], [1 : 0 : 1], [1 : 0 : 3], [1 : 1 : 0],
-         [1 : 1 : 1], [1 : 1 : 2], [1 : 1 : 3], [1 : 1 : 4],
-         [1 : 1 : 5], [1 : 1 : 6], [1 : 3 : 0], [1 : 3 : 1],
-         [1 : 3 : 2], [1 : 3 : 3], [1 : 3 : 4], [1 : 3 : 5],
-         [1 : 3 : 6]]
+        [[0 : 0 : 1],
+         [0 : 1 : 0],
+         [0 : 1 : 1],
+         [0 : 1 : 5],
+         [1 : 0 : 0],
+         [1 : 0 : 1],
+         [1 : 0 : 5],
+         [1 : 1 : 0],
+         [1 : 1 : 1],
+         [1 : 1 : 2],
+         [1 : 1 : 3],
+         [1 : 1 : 4],
+         [1 : 1 : 5],
+         [1 : 1 : 6],
+         [1 : 3 : 1],
+         [1 : 3 : 2],
+         [1 : 3 : 3],
+         [1 : 3 : 4],
+         [1 : 3 : 5],
+         [1 : 3 : 6],
+         [1 : 5 : 0]]
 
     As for a non-compact example, the blow-up of the plane is the line
     bundle $O_{\mathbf{P}^1}(-1)$. Its point set is the Cartesian
@@ -641,7 +656,7 @@ class SchemeHomset_points_subscheme_toric_field(SchemeHomset_points_toric_base):
             sage: P2.<x,y,z> = toric_varieties.P2(base_ring=GF(5))
             sage: cubic = P2.subscheme([x^3 + y^3 + z^3])
             sage: list(cubic.point_set())
-            [[0 : 1 : 4], [1 : 0 : 4], [1 : 4 : 0], [1 : 2 : 1], [1 : 1 : 2], [1 : 3 : 3]]
+            [[0 : 1 : 4], [1 : 0 : 4], [1 : 4 : 0], [1 : 1 : 2], [1 : 2 : 1], [1 : 3 : 3]]
             sage: cubic.point_set().cardinality()
             6
         """
@@ -661,7 +676,7 @@ class SchemeHomset_points_subscheme_toric_field(SchemeHomset_points_toric_base):
             sage: P2.<x,y,z> = toric_varieties.P2(base_ring=GF(5))
             sage: cubic = P2.subscheme([x^3 + y^3 + z^3])
             sage: list(cubic.point_set())
-            [[0 : 1 : 4], [1 : 0 : 4], [1 : 4 : 0], [1 : 2 : 1], [1 : 1 : 2], [1 : 3 : 3]]
+            [[0 : 1 : 4], [1 : 0 : 4], [1 : 4 : 0], [1 : 1 : 2], [1 : 2 : 1], [1 : 3 : 3]]
             sage: cubic.point_set().cardinality()
             6
         """
diff --git a/src/sage/schemes/toric/points.py b/src/sage/schemes/toric/points.py
index 31e7769ede..361a010d2b 100644
--- a/src/sage/schemes/toric/points.py
+++ b/src/sage/schemes/toric/points.py
@@ -537,7 +537,7 @@ class FiniteFieldPointEnumerator(NaiveFinitePointEnumerator):
             sage: enum._Chow_group_torsion()
             ((1, 2, 4), (1, 4, 2))
             sage: enum._Chow_group_torsion_generators()
-            ((1, 2, 4),)
+            ((1, 4, 2),)
         """
         if self.fan.is_smooth():
             return tuple()
@@ -673,7 +673,7 @@ class FiniteFieldPointEnumerator(NaiveFinitePointEnumerator):
             sage: list(cokernel)
             [(0), (1)]
             sage: [p.lift() for p in cokernel]
-            [(0, 0), (0, 1)]
+            [(0, 0), (0, -1)]
         """
         from sage.matrix.constructor import matrix, block_matrix, identity_matrix
         from sage.rings.all import ZZ
@@ -955,9 +955,9 @@ class FiniteFieldSubschemePointEnumerator(NaiveSubschemePointEnumerator):
             sage: ffe.homogeneous_coordinates([0], nonzero_coordinates, cokernel)
             (1, 1, 0)
             sage: ffe.homogeneous_coordinates([1], nonzero_coordinates, cokernel)
-            (1, 3, 0)
+            (1, 5, 0)
             sage: ffe.homogeneous_coordinates([2], nonzero_coordinates, cokernel)
-            (1, 2, 0)
+            (1, 4, 0)
         """
         z = [self.ambient.ring.zero()] * len(self.ambient.rays())
         z_nonzero = self.ambient.exp(
@@ -986,7 +986,7 @@ class FiniteFieldSubschemePointEnumerator(NaiveSubschemePointEnumerator):
             sage: point_set = X.point_set()
             sage: ffe = point_set._enumerator()
             sage: list(ffe)   # indirect doctest
-            [(1, 4, 3), (1, 1, 6), (1, 2, 5)]
+            [(1, 1, 6), (1, 2, 5), (1, 4, 3)]
         """
         for cone, nonzero_coordinates, cokernel in self.ambient.cone_points_iter():
             R = PolynomialRing(self.ambient.ring, cokernel.ngens(), 't')
@@ -1011,11 +1011,11 @@ class FiniteFieldSubschemePointEnumerator(NaiveSubschemePointEnumerator):
             sage: Y = X.subscheme(u^3 + v^3 + w^3 + u*v*w)
             sage: point_set = Y.point_set()
             sage: list(point_set)
-            [[0 : 1 : 3],
-             [1 : 0 : 3],
-             [1 : 3 : 0],
-             [1 : 1 : 6],
+            [[0 : 1 : 5],
+             [1 : 0 : 5],
+             [1 : 5 : 0],
              [1 : 1 : 4],
+             [1 : 1 : 6],
              [1 : 3 : 2],
              [1 : 3 : 5]]
             sage: ffe = point_set._enumerator()
diff --git a/src/sage/structure/factorization.py b/src/sage/structure/factorization.py
index 1d32db0842..7636f1a9ba 100644
--- a/src/sage/structure/factorization.py
+++ b/src/sage/structure/factorization.py
@@ -133,17 +133,17 @@ Factorizations can involve fairly abstract mathematical objects::
     sage: K.<a> = NumberField(x^2 + 3); K
     Number Field in a with defining polynomial x^2 + 3
     sage: f = K.factor(15); f
-    (Fractional ideal (-a))^2 * (Fractional ideal (5))
+    (Fractional ideal (1/2*a + 3/2))^2 * (Fractional ideal (5))
     sage: f.universe()
     Monoid of ideals of Number Field in a with defining polynomial x^2 + 3
     sage: f.unit()
     Fractional ideal (1)
     sage: g=K.factor(9); g
-    (Fractional ideal (-a))^4
+    (Fractional ideal (1/2*a + 3/2))^4
     sage: f.lcm(g)
-    (Fractional ideal (-a))^4 * (Fractional ideal (5))
+    (Fractional ideal (1/2*a + 3/2))^4 * (Fractional ideal (5))
     sage: f.gcd(g)
-    (Fractional ideal (-a))^2
+    (Fractional ideal (1/2*a + 3/2))^2
     sage: f.is_integral()
     True
 
diff --git a/src/sage/tests/books/computational-mathematics-with-sagemath/linalg_doctest.py b/src/sage/tests/books/computational-mathematics-with-sagemath/linalg_doctest.py
index d67a33e7c4..7cba0ef6bb 100644
--- a/src/sage/tests/books/computational-mathematics-with-sagemath/linalg_doctest.py
+++ b/src/sage/tests/books/computational-mathematics-with-sagemath/linalg_doctest.py
@@ -232,13 +232,13 @@ Sage example in ./linalg.tex, line 1640::
   sage: A = matrix(ZZ, 4, 5,\
   ....:            [-1,-1,-1,-2,-2,-2,1,1,-1,2,2,2,2,2,-1,2,2,2,2,2])
   sage: S,U,V = A.smith_form(); S,U,V
-  (
-                              [ 0 -2 -1 -5  0]
-  [1 0 0 0 0]  [ 1  0  0  0]  [ 1  0  1 -1 -1]
-  [0 1 0 0 0]  [ 0  0  1  0]  [ 0  0  0  0  1]
-  [0 0 3 0 0]  [-2  1  0  0]  [-1  2  0  5  0]
-  [0 0 0 6 0], [ 0  0 -2 -1], [ 0 -1  0 -2  0]
-  )
+    (
+                                [ 0 -2  0 -5  0]
+    [1 0 0 0 0]  [ 1  0  0  0]  [ 1  0  0 -1 -1]
+    [0 1 0 0 0]  [-2  1  1  0]  [ 0  0  0  0  1]
+    [0 0 3 0 0]  [-4  2  3  0]  [-1  1  1  5  0]
+    [0 0 0 6 0], [ 4 -2 -2 -1], [ 0  0 -1 -2  0]
+    )
 
 Sage example in ./linalg.tex, line 1674::
 
diff --git a/src/sage/tests/books/judson-abstract-algebra/galois-sage.py b/src/sage/tests/books/judson-abstract-algebra/galois-sage.py
index 6c25aa5dc5..8bdceb9424 100644
--- a/src/sage/tests/books/judson-abstract-algebra/galois-sage.py
+++ b/src/sage/tests/books/judson-abstract-algebra/galois-sage.py
@@ -385,69 +385,49 @@ r"""
 
     sage: L.subfields()
     [
-    (Number Field in c0 with defining polynomial x,
-     Ring morphism:
-       From: Number Field in c0 with defining polynomial x
-       To:   Number Field in c with defining polynomial x^8 + 28*x^4 + 2500
-       Defn: 0 |--> 0,
-     None),
-    (Number Field in c1 with defining polynomial x^2 + 112*x + 40000,
-     Ring morphism:
-       From: Number Field in c1 with defining polynomial x^2 + 112*x + 40000
-       To:   Number Field in c with defining polynomial x^8 + 28*x^4 + 2500
-       Defn: c1 |--> 4*c^4,
-     None),
-    (Number Field in c2 with defining polynomial x^2 + 512,
-     Ring morphism:
-       From: Number Field in c2 with defining polynomial x^2 + 512
-       To:   Number Field in c with defining polynomial x^8 + 28*x^4 + 2500
-       Defn: c2 |--> 1/25*c^6 + 78/25*c^2,
-     None),
-    (Number Field in c3 with defining polynomial x^2 - 288,
-     Ring morphism:
-       From: Number Field in c3 with defining polynomial x^2 - 288
-       To:   Number Field in c with defining polynomial x^8 + 28*x^4 + 2500
-       Defn: c3 |--> -1/25*c^6 + 22/25*c^2,
-     None),
-    (Number Field in c4 with defining polynomial x^4 + 112*x^2 + 40000,
-     Ring morphism:
-       From: Number Field in c4 with defining polynomial x^4 + 112*x^2 + 40000
-       To:   Number Field in c with defining polynomial x^8 + 28*x^4 + 2500
-       Defn: c4 |--> 2*c^2,
-     None),
-    (Number Field in c5 with defining polynomial x^4 + 648,
-     Ring morphism:
-       From: Number Field in c5 with defining polynomial x^4 + 648
-       To:   Number Field in c with defining polynomial x^8 + 28*x^4 + 2500
-       Defn: c5 |--> 1/80*c^5 + 79/40*c,
-     None),
-    (Number Field in c6 with defining polynomial x^4 + 8,
-    Ring morphism:
-      From: Number Field in c6 with defining polynomial x^4 + 8
+    (Number Field in c0 with defining polynomial x, Ring morphism:
+      From: Number Field in c0 with defining polynomial x
       To:   Number Field in c with defining polynomial x^8 + 28*x^4 + 2500
-      Defn: c6 |--> -1/80*c^5 + 1/40*c,
-      None),
-    (Number Field in c7 with defining polynomial x^4 - 512,
-     Ring morphism:
-       From: Number Field in c7 with defining polynomial x^4 - 512
-       To:   Number Field in c with defining polynomial x^8 + 28*x^4 + 2500
-       Defn: c7 |--> -1/60*c^5 + 41/30*c,
-     None),
-    (Number Field in c8 with defining polynomial x^4 - 32,
-     Ring morphism:
-       From: Number Field in c8 with defining polynomial x^4 - 32
-       To:   Number Field in c with defining polynomial x^8 + 28*x^4 + 2500
-       Defn: c8 |--> 1/60*c^5 + 19/30*c,
-     None),
-    (Number Field in c9 with defining polynomial x^8 + 28*x^4 + 2500,
-     Ring morphism:
-       From: Number Field in c9 with defining polynomial x^8 + 28*x^4 + 2500
-       To:   Number Field in c with defining polynomial x^8 + 28*x^4 + 2500
-       Defn: c9 |--> c,
-     Ring morphism:
-       From: Number Field in c with defining polynomial x^8 + 28*x^4 + 2500
-       To:   Number Field in c9 with defining polynomial x^8 + 28*x^4 + 2500
-       Defn: c |--> c9)
+      Defn: 0 |--> 0, None),
+    (Number Field in c1 with defining polynomial x^2 + 112*x + 40000, Ring morphism:
+      From: Number Field in c1 with defining polynomial x^2 + 112*x + 40000
+      To:   Number Field in c with defining polynomial x^8 + 28*x^4 + 2500
+      Defn: c1 |--> 4*c^4, None),
+    (Number Field in c2 with defining polynomial x^2 + 512, Ring morphism:
+      From: Number Field in c2 with defining polynomial x^2 + 512
+      To:   Number Field in c with defining polynomial x^8 + 28*x^4 + 2500
+      Defn: c2 |--> 1/25*c^6 + 78/25*c^2, None),
+    (Number Field in c3 with defining polynomial x^2 - 288, Ring morphism:
+      From: Number Field in c3 with defining polynomial x^2 - 288
+      To:   Number Field in c with defining polynomial x^8 + 28*x^4 + 2500
+      Defn: c3 |--> -1/25*c^6 + 22/25*c^2, None),
+    (Number Field in c4 with defining polynomial x^4 + 112*x^2 + 40000, Ring morphism:
+      From: Number Field in c4 with defining polynomial x^4 + 112*x^2 + 40000
+      To:   Number Field in c with defining polynomial x^8 + 28*x^4 + 2500
+      Defn: c4 |--> 2*c^2, None),
+    (Number Field in c5 with defining polynomial x^4 + 8, Ring morphism:
+      From: Number Field in c5 with defining polynomial x^4 + 8
+      To:   Number Field in c with defining polynomial x^8 + 28*x^4 + 2500
+      Defn: c5 |--> -1/80*c^5 + 1/40*c, None),
+    (Number Field in c6 with defining polynomial x^4 + 648, Ring morphism:
+      From: Number Field in c6 with defining polynomial x^4 + 648
+      To:   Number Field in c with defining polynomial x^8 + 28*x^4 + 2500
+      Defn: c6 |--> 1/80*c^5 + 79/40*c, None),
+    (Number Field in c7 with defining polynomial x^4 - 512, Ring morphism:
+      From: Number Field in c7 with defining polynomial x^4 - 512
+      To:   Number Field in c with defining polynomial x^8 + 28*x^4 + 2500
+      Defn: c7 |--> -1/60*c^5 + 41/30*c, None),
+    (Number Field in c8 with defining polynomial x^4 - 32, Ring morphism:
+      From: Number Field in c8 with defining polynomial x^4 - 32
+      To:   Number Field in c with defining polynomial x^8 + 28*x^4 + 2500
+      Defn: c8 |--> 1/60*c^5 + 19/30*c, None),
+    (Number Field in c9 with defining polynomial x^8 + 28*x^4 + 2500, Ring morphism:
+      From: Number Field in c9 with defining polynomial x^8 + 28*x^4 + 2500
+      To:   Number Field in c with defining polynomial x^8 + 28*x^4 + 2500
+      Defn: c9 |--> c, Ring morphism:
+      From: Number Field in c with defining polynomial x^8 + 28*x^4 + 2500
+      To:   Number Field in c9 with defining polynomial x^8 + 28*x^4 + 2500
+      Defn: c |--> c9)
     ]
 
 ~~~~~~~~~~~~~~~~~~~~~~ ::
diff --git a/src/sage/tests/parigp.py b/src/sage/tests/parigp.py
index c118b6eb37..f567ade058 100644
--- a/src/sage/tests/parigp.py
+++ b/src/sage/tests/parigp.py
@@ -2,15 +2,6 @@ r"""
 This file is meant to catch errors in the PARI/GP package which are not
 caught by any other tests.
 
-Check that :trac:`9876` has been fixed, this test comes from PARI's
-self-test :pari:`rnfkummer` but was modified such that the answer is
-canonical::
-
-    sage: pari('addprimes([31438243, 238576291, 18775387483, 24217212463267, 1427657500359111961, 135564809928627323997297867381959])')
-    [31438243, 238576291, 18775387483, 24217212463267, 1427657500359111961, 135564809928627323997297867381959]
-    sage: pari('K = bnfinit(y^4-52*y^2+26,1); pol = rnfkummer(bnrinit(K,3,1),Mat(5)); L = rnfinit(K, pol); polredabs(polredbest(L.polabs))')  # long time
-    x^20 - 112*x^18 + 5108*x^16 - 123460*x^14 + 1724337*x^12 - 14266996*x^10 + 69192270*x^8 - 188583712*x^6 + 260329852*x^4 - 141461008*x^2 + 19860776
-
 Check that :trac:`10195` (PARI bug 1153) has been fixed::
 
     sage: print(gp.eval("mathnf([0,0,0,0,0,0,0,0,0,13;0,0,0,0,0,0,0,0,23,6;0,0,0,0,0,0,0,23,-4,-7;0,0,0,0,0,0,17,-3,5,-5;0,0,0,0,0,56,16,-16,-15,-17;0,0,0,0,57,24,-16,-25,2,-21;0,0,0,114,9,56,51,-52,25,-55;0,0,113,-31,-11,24,0,28,34,-16;0,50,3,2,16,-6,-2,7,-19,-21;118,43,51,23,37,-52,18,38,51,28],0)"))