Primitive element is slow for RootOf

[oscarbenjamin]

This is needed e.g. to construct a domain for the Jordan normal form of a matrix.

I'm wondering if there is an algorithm that can do this more directly. Specifically given an irreducible polynomial is there an efficient way to construct a primitive element for the splitting field (generated by the roots)?

We can compute a primitive element for the splitting field of an irreducible polynomial as:

from sympy import *
from sympy.abc import x

rs = [rootof(2*x**3 + x + 2, n) for n in range(3)]

# Takes around 1 minute:
p, inc, *ex = primitive_element(rs, ex=True)

print(p)
print(inc)
print(ex)

However this is very slow. We can make it a bit faster with (see also #19732):

diff --git a/sympy/polys/numberfields.py b/sympy/polys/numberfields.py
index 79f60d3e89..1a7bb2f84d 100644
--- a/sympy/polys/numberfields.py
+++ b/sympy/polys/numberfields.py
@@ -58,25 +58,21 @@ def _choose_factor(factors, x, v, dom=QQ, prec=200, bound=5):
     t = QQ(1, 10)
 
     for n in range(bound**len(symbols)):
-        prec1 = 10
         n_temp = n
         for s in symbols:
             points[s] = n_temp % bound
             n_temp = n_temp // bound
 
-        while True:
-            candidates = []
-            eps = t**(prec1 // 2)
-            for f in factors:
-                if abs(f.as_expr().evalf(prec1, points)) < eps:
-                    candidates.append(f)
-            if candidates:
-                factors = candidates
-            if len(factors) == 1:
-                return factors[0]
-            if prec1 > prec:
-                break
-            prec1 *= 2
+        def nonzero(f):
+            n10 = abs(f.as_expr()).evalf(10, points)
+            if n10._prec > 1:
+                return n10 > 0
+
+        candidates = [f for f in factors if not nonzero(f)]
+        if candidates:
+            factors = candidates
+        if len(factors) == 1:
+            return factors[0]
 
     raise NotImplementedError("multiple candidates for the minimal polynomial of %s" % v)
 
diff --git a/sympy/polys/rootisolation.py b/sympy/polys/rootisolation.py
index c2b44ef0a8..e71f10bb72 100644
--- a/sympy/polys/rootisolation.py
+++ b/sympy/polys/rootisolation.py
@@ -151,7 +151,8 @@ def dup_step_refine_real_root(f, M, K, fast=False):
     if a == b and c == d:
         return f, (a, b, c, d)
 
-    A = dup_root_lower_bound(f, K)
+    # A = dup_root_lower_bound(f, K)
+    A = None
 
     if A is not None:
         A = K(int(A))