GaloisGrp.jl

module GaloisGrp

using Oscar, Random
import Base: ^, +, -, *, ==
import Oscar: Hecke, AbstractAlgebra, GAP, extension_field, isinteger
using Oscar: SLPolyRing, SLPoly, SLPolynomialRing, CycleType

import Oscar: pretty, LowercaseOff

export cauchy_ideal
export elementary_symmetric
export fixed_field
export galois_group
export galois_ideal
export galois_quotient
export power_sum
export slpoly_ring
export to_elementary_symmetric
export upper_bound
export valuation_of_roots

import Hecke: orbit, fixed_field, extension_field

function __init__()
  GAP.Packages.load("ferret"; install=true)

  Hecke.add_verbosity_scope(:GaloisGroup)
  Hecke.add_verbosity_scope(:GaloisInvariant)
  Hecke.add_assertion_scope(:GaloisInvariant)
end

"""
A poor mans version of a (range of) tropical rings.

Actually, not even a ring. Defined in terms of operations
 - `mul` used for multiplication
 - `add` for addition
 - `pow` for powering with Int exponent
 - `map` to create elements from (large) integers/ ZZRingElem
 - `name` is only used for printing.
"""
struct BoundRing{T}  <: AbstractAlgebra.Ring
  mul#::(T,T) -> T
  add#::(T,T) -> T
  pow#::(T, Int) -> T
  map#:: R -> T
  name::String
end

function Base.show(io::IO, b::BoundRing{T}) where {T}
  print(io, "$(b.name) for type $T")
end

struct BoundRingElem{T} <: AbstractAlgebra.RingElem
  val::T
  p::BoundRing # the parent
end

function Base.show(io::IO, b::BoundRingElem)
  print(io, "(x <= $(b.val))")
end

function check_parent(a::BoundRingElem, b::BoundRingElem)
  parent(a) == parent(b) || error("Elements must have same parent")
  return true
end

Base.:(==)(a::BoundRingElem, b::BoundRingElem) = check_parent(a, b) && a.val == b.val
function +(a::BoundRingElem, b::BoundRingElem) 
  check_parent(a, b) 
  c = BoundRingElem(a.p.add(a.val, b.val), a.p)
#  @show a, "+", b, "=", c
  return c
end
-(a::BoundRingElem, b::BoundRingElem) = check_parent(a, b) && BoundRingElem(a.p.add(a.val, b.val), a.p)
function *(a::BoundRingElem, b::BoundRingElem) 
  check_parent(a, b) 
  c = BoundRingElem(a.p.mul(a.val, b.val), a.p)
#  @show a, "*", b, ":=", c
  return c
end

function *(a::ZZRingElem, b::BoundRingElem) 
  c = BoundRingElem(b.p.mul(b.p.map(a), b.val), b.p)
#  @show a, ":*", b, ":=", c
  return c
end
function ^(a::BoundRingElem, b::Int) 
  c = BoundRingElem(a.p.pow(a.val, b), a.p)
#  @show a, ":^", b, ":=", c
  return c
end
-(a::BoundRingElem) = a

Oscar.parent(a::BoundRingElem) = a.p
value(a::BoundRingElem) = a.val
Base.isless(a::BoundRingElem, b::BoundRingElem) = check_parent(a, b) && isless(value(a), value(b))

Oscar.parent_type(::Type{BoundRingElem{T}}) where T = BoundRing{T}
Oscar.elem_type(::Type{BoundRing{T}}) where T = BoundRingElem{T}

(R::BoundRing{T})(a::ZZRingElem) where T = BoundRingElem{T}(R.map(abs(a)), R)
(R::BoundRing{T})(a::Integer) where T = BoundRingElem{T}(ZZRingElem(a), R)
(R::BoundRing{T})() where T = BoundRingElem{T}(ZZRingElem(0), R)
(R::BoundRing{T})(a::T) where T = BoundRingElem{T}(a, R)
(R::BoundRing{T})(a::BoundRingElem{T}) where T = a
Oscar.one(R::BoundRing) = R(1)
Oscar.zero(R::BoundRing) = R(0)


"""
Intended to be used to get upper bounds under evaluation in function fields (non 
  Archimedean valuaton, degree)

Operations are: 
 - `a+b := max(a, b)`
 - `ab  := a+b`
"""
function max_ring()
  return BoundRing{ZZRingElem}( (x,y) -> x+y, (x,y) -> max(x, y), (x,y) -> y*x, x->x, "max-ring")
end

"""
Normal ring
"""
function add_ring(;type::Type=ZZRingElem)
  return BoundRing{type}( (x,y) -> x*y, (x,y) -> x+y, (x,y) -> x^y, x->abs(x), "add-ring")
end

#roots rt are power series sum a_n x^n
#we have |a_n| <= r^-n B (n+1)^k for B = x[1], k = x[2]
#and deg(rt) <= x[3] (infinite valuation)
function qt_ring()
  return BoundRing{Tuple{ZZRingElem, Int, QQFieldElem}}( (x,y) -> (x[1]*y[1], x[2]+y[2]+1, x[3]+y[3]),
                                            (x,y) -> (x[1]+y[1], max(x[2], y[2]), max(x[3], y[3])),
                                            (x,y) -> (x[1]^y, y*x[2]+y-1, y*x[3]),
                                                x -> _coerce_qt(x), "qt-ring")
end

_coerce_qt(x::ZZRingElem) = (abs(x), 0, QQFieldElem(0))
_coerce_qt(x::Integer) = (ZZRingElem(x), 0, QQFieldElem(0))
function _coerce_qt(x::ZZPolyRingElem)
  if iszero(x)
    return (ZZRingElem(0), 0, QQFieldElem(0))
  end
  return (maximum(abs, coefficients(x))*(degree(x)+1), 0, QQFieldElem(0))
end

(R::BoundRing{Tuple{ZZRingElem, Int, QQFieldElem}})(a::Tuple{ZZRingElem, Int, QQFieldElem}) = BoundRingElem(a, R)
(R::BoundRing{Tuple{ZZRingElem, Int, QQFieldElem}})(a::Integer) = BoundRingElem(R.map(a), R)
(R::BoundRing{Tuple{ZZRingElem, Int, QQFieldElem}})(a::ZZRingElem) = BoundRingElem(R.map(a), R)
(R::BoundRing{Tuple{ZZRingElem, Int, QQFieldElem}})(a::ZZPolyRingElem) = BoundRingElem(R.map(a), R)

"""
An slpoly evaluated at `cost_ring` elements `0` will count the number
of multiplications involved. A measure of the cost of evaluation at more
interesting scalars.

Operations:
 - `xy := x+y+1`
 - `x+y := x+y`
 - `x^y := x+2*log_2(y)`
 - all constants are mapped to `0`
"""
function cost_ring()
  return BoundRing{ZZRingElem}( (x,y) -> x+y+1, (x,y) -> x+y, (x,y) -> x+2*nbits(y), x->0, "cost-ring")
end

"""
An slpoly evaluated at `degree_ring` elements `1` will bound the total degree
from above.

Operations:
 - `xy := x+y`
 - `x+y := max(x,y)`
 - `x^y := yx`
 - all constants are mapped to `0`
"""
function degree_ring()
  return BoundRing{ZZRingElem}( (x,y) -> x+y, (x,y) -> max(x, y), (x,y) -> y*x, x->0, "degree-ring")
end

@doc raw"""
    cost(I::SLPoly)

Counts the number of multiplications to evaluate `I`, optionally
a Tschirnhaus transformation (`ZZPolyRingElem`) can be passed in as well.
"""
function cost(I::SLPoly)
  n = ngens(parent(I))
  C = cost_ring()
  return value(evaluate(I, [C(0) for i = 1:n]))
end
function cost(I::SLPoly, ts::ZZPolyRingElem)
  n = ngens(parent(I))
  C = cost_ring()
  return value(evaluate(I, [C(0) for i = 1:n]))+n*degree(ts)
end

@doc raw"""
    total_degree(I::SLPoly)

Determines an upper bound for the total degree of `I`.
"""
function total_degree(I::SLPoly)
  n = ngens(parent(I))
  C = degree_ring()
  return value(evaluate(I, [C(1) for i = 1:n]))
end

Oscar.mul!(a::BoundRingElem, b::BoundRingElem, c::BoundRingElem) = b*c
Oscar.addeq!(a::BoundRingElem, b::BoundRingElem) = a+b

#my 1st invariant!!!
@doc raw"""
    sqrt_disc(a::Vector)

The product of differences ``a[i] - a[j]`` for all indices ``i<j``.    
"""
function sqrt_disc(a::Vector)
  if length(a) == 1
    return one(parent(a[1]))
  end
  return prod([a[i] - a[j] for i = 1:length(a)-1 for j = i+1:length(a)])
end

@doc raw"""
    elementary_symmetric(g::Vector, i::Int)

Evaluates the `i`-th elementary symmetric polynomial at the values in `g`.    
The `i`-th elementary symmetric polynomial is the sum over all
products of `i` distinct variables.
"""
function elementary_symmetric(g::Vector, i::Int)
  return sum(prod(g[i] for i = s) for s = Hecke.subsets(Set(1:length(g)), i))
end

@doc raw"""
    power_sum(g::Vector, i::Int)

Evaluates the `i`-th power sums at the values in `g`, ie. the sum
of the `i`-th power of the values.
"""
function power_sum(g::Vector, i::Int)
  return sum(a^i for a = g)
end

@doc raw"""
    discriminant(g::Vector)

Compute the product of all differences of distinct elements in the array.    
"""
function Oscar.discriminant(g::Vector{<:RingElem})
  return prod(a-b for a = g for b = g if a!=b)
end


function slpoly_ring(R::AbstractAlgebra.Ring, n::Int; cached::Bool = false)
  return SLPolynomialRing(R, [ Symbol("x_$i") for i=1:n], cached = cached)
end

function slpoly_ring(R::AbstractAlgebra.Ring, p::Pair{Symbol, <:AbstractVector{Int}}...; cached::Bool = false)
  return SLPolynomialRing(R, p..., cached = cached)
end

function (R::SLPolyRing)(a::SLPoly)
  parent(a) == R && return a
  error("wrong parent")
end

@doc raw"""
    roots_upper_bound(f::ZZPolyRingElem) -> ZZRingElem

An upper upper_bound for the absolute value of the complex roots of the input.    
Uses the Cauchy bound.
"""
function Nemo.roots_upper_bound(f::ZZPolyRingElem)
  a = coeff(f, degree(f))
  b = ceil(ZZRingElem, abs(coeff(f, degree(f)-1)//a))
  for i=0:degree(f)-2
    b = max(b, iroot(ceil(ZZRingElem, abs(coeff(f, i)//a)), degree(f)-i)+1)
  end
  return 2*b
  return max(ZZRingElem(1), maximum([ceil(ZZRingElem, abs(coeff(f, i)//a)) for i=0:degree(f)]))
end
function Nemo.roots_upper_bound(f::QQPolyRingElem)
  a = coeff(f, degree(f))
  return max(ZZRingElem(1), maximum([ceil(ZZRingElem, abs(coeff(f, i)//a)) for i=0:degree(f)]))
end

#roots are sums of m distinct roots of f
#from https://doi.org/10.2307/2153295
#Symmetric Functions, m-Sets, and Galois Groups
#by David Casperson and John McKay
@doc raw"""
    msum_poly(f::PolyRingElem, m::Int)

Compute the polynomial with roots sums of `m` roots of `f` using
resultants.
"""
function msum_poly(f::PolyRingElem, m::Int)
  f = divexact(f, leading_coefficient(f))
  N = binomial(degree(f), m)
  p = Hecke.polynomial_to_power_sums(f, N)
  p = vcat([degree(f)*one(base_ring(f))], p)
  S, a = power_series_ring(base_ring(f), N+1, "a")
  Hfs = S([p[i]//factorial(ZZRingElem(i-1)) for i=1:length(p)], N+1, N+1, 0)
  H = [S(1), Hfs]
  for i=2:m
    push!(H, 1//i*sum((-1)^(h+1)*Hfs(h*a)*H[i-h+1] for h=1:i))
  end
  p = [coeff(H[end], i)*factorial(ZZRingElem(i)) for i=0:N]
  return Hecke.power_sums_to_polynomial(p[2:end])
end

@doc raw"""
A `GaloisCtx`, is the context object used for the computation of Galois
groups of (univariate) polynomials. It contains
 - the polynomial
 - an object that can compute the roots in a fixed order up to a given 
   precision. Currently, this is a q-adic field, that is an unramified
   extension of the p-adic numbers.
 - a upper bound on the _size_ of the roots, currently an upper bound on the
   complex absolute value.
 - at the end, the Galois group  

This is constructed implicitly while computing a Galois group and returned
together with the group.

Not type stable, not sure what to do about it:
Depends on type of
 - `f`
However, `f` is hardly ever used. 
"""
mutable struct GaloisCtx{T}
  f::PolyRingElem
  C::T # a suitable root context
  B::BoundRingElem # a "bound" on the roots, might be "anything"
  G::PermGroup
  rt_num::Dict{Int, Int}
  chn::Vector{Tuple{PermGroup, SLPoly, ZZPolyRingElem, Vector{PermGroupElem}}}
  start::Tuple{Int, Vector{Vector{Vector{Int}}}} # data for the starting group:
    # if start[1] == 1: start[2] is a list of the block systems used
    #             == 2  start[2] is a list of the orbits of the msum-poly, ie. a list of a list of pairs
    #                            where the pairs are Vector{Int} of length 2
  data::Any #whatever else is needed in special cases
  #= the descent chain, recording
   - the group
   - the invariant
   - the tschirnhaus transformation
   - the cosets used
   should probably also record if the step was proven or not
    the starting group  and the block systems used to get them
  =#
  prime::Any #=can be
   - ZZRingElem/ Int: prime number, used over Q
   - AbsSimpleNumFieldOrderIdeal : prime ideal , used over NfAbs
   - (Int, Int): evaluation point, prime number used over Q(t)
   =#

  function GaloisCtx(f::ZZPolyRingElem, ::AcbField)
    r = new{ComplexRootCtx}()
    r.f = f
    r.C = ComplexRootCtx(f)
    r.B = add_ring()(leading_coefficient(f)*roots_upper_bound(f))
    r.chn = Tuple{PermGroup, SLPoly, ZZPolyRingElem, Vector{PermGroupElem}}[]
    return r
  end

  function GaloisCtx(f::QQPolyRingElem, x::AcbField)
    return GaloisCtx(numerator(f), x)
  end

  function GaloisCtx(f::ZZPolyRingElem, p::Int)
    r = new{Hecke.qAdicRootCtx}()
    r.prime = p
    r.f = f
    r.C = Hecke.qAdicRootCtx(f, p, splitting_field = true)
    r.B = add_ring()(leading_coefficient(f)*roots_upper_bound(f))
    r.chn = Tuple{PermGroup, SLPoly, ZZPolyRingElem, Vector{PermGroupElem}}[]
    return r
  end

  function GaloisCtx(f::ZZPolyRingElem, field::Union{Nothing, AbsSimpleNumField})
    r = new{SymbolicRootCtx}()
    r.f = f
    r.C = SymbolicRootCtx(f, field)
    r.B = add_ring()(leading_coefficient(f)*roots_upper_bound(f))
    r.chn = Tuple{PermGroup, SLPoly, ZZPolyRingElem, Vector{PermGroupElem}}[]
    return r
  end

  function GaloisCtx(f::QQPolyRingElem, p::Int)
    d = mapreduce(denominator, lcm, coefficients(f))
    return GaloisCtx(Hecke.Globals.Zx(d*f), p)
  end
  #=
  Roots in F_q[[t]] for q = p^d
    - needs to be tweaked to do Q_q[[t]
    - need q-lifting as well as t-lifting
    - possibly also mul_ks for Q_q[[t]] case
    - needs merging in Hecke
  =#

  function GaloisCtx(f::ZZMPolyRingElem, shft::Int, p::Int, d::Int)
    f = evaluate(f, [gen(parent(f), 1), gen(parent(f), 2)+shft])
    #f(x, T+t), the roots are power series in T over qAdic(p, d)
    #so basically for f in Qq<<(T+t)>> 
    @assert ngens(parent(f)) == 2
    Qq, _ = QadicField(p, d, 10)
    F, mF = residue_field(Qq)
    H = Hecke.MPolyFact.HenselCtxFqRelSeries(f, F)
    SQq, _ = power_series_ring(Qq, 2, "s", cached = false)
    SQqt, _ = polynomial_ring(SQq, "t", cached = false)
    mc(f) = map_coefficients(x->map_coefficients(y->setprecision(preimage(mF, y), 1), x, parent = SQq), f, parent = SQqt)
    HQ = Hecke.MPolyFact.HenselCtxFqRelSeries(H.f, map(mc, H.lf), map(mc, H.cf), H.n)
    r = new{Hecke.MPolyFact.HenselCtxFqRelSeries{AbstractAlgebra.Generic.RelSeries{QadicFieldElem}}}()
    r.prime = (shft, p)
    Qt, t = rational_function_field(QQ, "t", cached = false)
    Qts, s = polynomial_ring(Qt, "s", cached = false)
    r.f = evaluate(f, [s, Qts(t)])
    r.C = HQ
    r.chn = Tuple{PermGroup, SLPoly, ZZPolyRingElem, Vector{PermGroupElem}}[]
    
    vl = roots_upper_bound(f)
    r.B = qt_ring()(vl[1])
    r.data = Any[vl[2], shft, false] #false: not simulating C
    @assert typeof(r.data[3]) == Bool
    return r
  end
  function GaloisCtx(T::Type)
    r = new{T}()
    r.chn = Tuple{PermGroup, SLPoly, ZZPolyRingElem, Vector{PermGroupElem}}[]
    return r
  end
end

function Oscar.prime(C::GaloisCtx{Hecke.MPolyFact.HenselCtxFqRelSeries{Generic.RelSeries{QadicFieldElem}}})
  return prime(base_ring(base_ring(C.C.lf[1])))
end

function bound_to_precision(G::GaloisCtx{T}, B::BoundRingElem{Tuple{ZZRingElem, Int, QQFieldElem}}, extra=(0, 0)) where {T}
  if isa(extra, Int)
    extra = (extra, min(2, div(extra, 3)))
  end
  C, k, d = B.val
  r = G.data[1]
  #so power series prec need to be floor(Int, d)
  n = floor(ZZRingElem, d+1)
  #padic: we ne |a_i| for i=0:n and |a_i| <= C (i+1)^k/r^i
  #and then log_p()
  #according to the Qt file, a_i is maximal around k/log(r) -1
  if G.data[3] #we're simulating CC, so we don't care about the p-adics (too much)
    b = max(C, ZZRingElem(10)^10)
  elseif isone(r)
    b = C*(n+1)^k
  else
    c = max(1, floor(Int, k/log(r)-1))
    if n<c
      b = C*(n+1)^k//r^n
    else
      b = max(C*(c+1)^k//r^c, C*(c+2)^k//r^(c+1))
    end
  end
  b = max(b, ZZRingElem(1))
  return (clog(floor(ZZRingElem, b), prime(G))+extra[1], Int(n)+extra[2])
end

function bound_to_precision(G::GaloisCtx{T}, B::BoundRingElem{ZZRingElem}, extra::Int = 5) where {T}
  return clog(B.val, G.C.p)+extra
end

mutable struct ComplexRootCtx
  f::ZZPolyRingElem
  pr::Int
  rt::Vector{AcbFieldElem}
  function ComplexRootCtx(f::ZZPolyRingElem)
    @assert is_monic(f)
    rt = roots(AcbField(20), f)
    return new(f, 20, rt)
  end
  function ComplexRootCtx(f::QQPolyRingElem)
    return ComplexRootCtx(numerator(f))
  end
end

function Base.show(io::IO, GC::GaloisCtx{ComplexRootCtx})
  print(pretty(io), LowercaseOff(), "Galois context for $(GC.f) using complex roots")
end

function Hecke.roots(C::GaloisCtx{ComplexRootCtx}, pr::Int = 10; raw::Bool = false)
  if C.C.pr >= pr
    return C.C.rt
  end
  rt = roots(AcbField(pr), C.C.f)
  C.C.pr = pr
  n = length(rt)
  for i=1:n
    C.C.rt[i] = rt[argmin([abs(C.C.rt[i] - rt[x]) for x = 1:n])]
  end
  return C.C.rt
end

function isinteger(GC::GaloisCtx{ComplexRootCtx}, B::BoundRingElem, e)
  if abs(imag(e)) > 1e-10
    return false, ZZRingElem(0)
  end
  r = round(ZZRingElem, real(e))
  if abs(real(e)-r) > 1e-10
    return false, ZZRingElem(0)
  else
    return true, r
  end
end

function bound_to_precision(G::GaloisCtx{ComplexRootCtx}, B::BoundRingElem{ZZRingElem}, extra::Int = 5)
  return 2*clog(B.val, 2) + 10
end

function map_coeff(G::GaloisCtx{ComplexRootCtx}, a::QQFieldElem)
  return parent(G.C.rt[1])(a)
end

function Hecke.MPolyFact.block_system(a::Vector{AcbFieldElem}, eps = 1e-9)
  b = Dict{Int, Vector{Int}}()
  for i=1:length(a)
    cb = collect(keys(b))
    fl = findfirst(x->abs(a[i] - a[x]) < eps, cb)
    if fl === nothing
      b[i] = [i]
    else
      push!(b[cb[fl]], i)
    end
  end
  bs = sort(collect(values(b)), lt = (a,b) -> isless(a[1], b[1]))
  return bs
end

mutable struct SymbolicRootCtx
  f::ZZPolyRingElem
  rt::Vector{AbsSimpleNumFieldElem}
  function SymbolicRootCtx(f::ZZPolyRingElem, ::Nothing)
    @assert is_monic(f)
    _, rt = splitting_field(f, do_roots = true)
    return new(f, rt)
  end
  function SymbolicRootCtx(f::ZZPolyRingElem, field::AbsSimpleNumField)
    @assert is_monic(f)
    rt = roots(f, field)
    return new(f, rt)
  end

  function SymbolicRootCtx(f::QQPolyRingElem)
    return SymbolicRootCtx(numerator(f), nothing)
  end
end

function Base.show(io::IO, GC::GaloisCtx{SymbolicRootCtx})
  print(pretty(io), LowercaseOff(), "Galois context for $(GC.f) using symbolic roots")
end

function Hecke.roots(C::GaloisCtx{SymbolicRootCtx}, ::Int; raw::Bool = false)
  return C.C.rt
end

function isinteger(GC::GaloisCtx{SymbolicRootCtx}, B::BoundRingElem, e)
  if Oscar.is_integer(e)
    return true, ZZ(e)
  else
    return false, ZZRingElem(0)
  end
end

function bound_to_precision(G::GaloisCtx{SymbolicRootCtx}, B::BoundRingElem{ZZRingElem}, extra::Int = 5)
  return 1
end

function map_coeff(G::GaloisCtx{SymbolicRootCtx}, a::QQFieldElem)
  return parent(G.C.rt[1])(a)
end

function Nemo.roots_upper_bound(f::ZZMPolyRingElem, t::Int = 0)
  @assert nvars(parent(f)) == 2
  Qs, s = rational_function_field(FlintQQ, "t", cached = false)
  Qsx, x = polynomial_ring(Qs, cached = false)
  F = evaluate(f, [x, Qsx(s)])
  dis = numerator(discriminant(F))
  @assert !iszero(dis(t))
  rt = roots(AcbField(20), dis)
  r = Hecke.lower_bound(minimum([abs(x-t) for x = rt]), ZZRingElem)
  @assert r > 0
  ff = map_coefficients(abs, f)
  C = roots_upper_bound(Hecke.Globals.Zx(map(x->evaluate(x, ZZRingElem[r, 0]), coefficients(ff, 2))))
  C1 = maximum(map(x->evaluate(x, ZZRingElem[r, 0]), coefficients(ff, 2)))
  #the infinite valuation... need Newton
  vl = valuations_of_roots(F)
  return (C+1, 0, maximum(x[1] for x = vl)), r
end

function Base.show(io::IO, GC::GaloisCtx{Hecke.qAdicRootCtx})
  print(pretty(io), LowercaseOff(), "Galois context for $(GC.f) and prime $(GC.C.p)")
end
function Base.show(io::IO, GC::GaloisCtx{<:Hecke.MPolyFact.HenselCtxFqRelSeries})
  print(pretty(io), LowercaseOff(), "Galois context for $(GC.f)")
end


#TODO: change pr to be a "bound_ring_elem": in the Qt case this has to handle
#      both power series prec as well as q-adic...
@doc raw"""
    roots(G::GaloisCtx, pr::Int)

The roots of the polynomial used to define the Galois context in the fixed order
used in the algorithm. The roots are returned up to a precision of `pr`
p-adic digits, thus they are correct modulo ``p^{pr}``

For non-monic polynomials the roots are scaled by the leading coefficient.
If `raw` is set to true, the scaling is omitted.
The bound in the `GaloisCtx` is also adjusted.
"""
function Hecke.roots(G::GaloisCtx{Hecke.qAdicRootCtx}, pr::Int=5; raw::Bool = false)
  a = Hecke.roots(G.C, pr)::Vector{QadicFieldElem}
  b = Hecke.expand(a, all = true, flat = false, degs = Hecke.degrees(G.C.H))::Vector{QadicFieldElem}
  if isdefined(G, :rt_num)
    b = [b[G.rt_num[i]] for i=1:length(G.rt_num)]
  end
  if raw
    return b
  else
    return leading_coefficient(G.f) .* b
  end
end

function Hecke.setprecision(a::Generic.RelSeries, p::Int)
  b = parent(a)(a.coeffs, min(length(a.coeffs), p), p+valuation(a), valuation(a))
end

#TODO: does not really work: the polynomial is not truncated, only
#      the entry for precision is updated: 1 + ... + s^40 + O(s^2)
function Hecke.setprecision!(G::GaloisCtx{<:Hecke.MPolyFact.HenselCtxFqRelSeries}, pr::Tuple{Int, Int})
  c = coeff(G.C.lf[1], 0)
  @assert precision(c) >= pr[2]
  @assert precision(coeff(c, 0)) >= pr[1]
  for f = G.C.lf
    Hecke.set_precision!(f, pr[2])
  end
  for f = G.C.cf  
    Hecke.set_precision!(f, pr[2])
  end
end

function Hecke.roots(G::GaloisCtx{<:Hecke.MPolyFact.HenselCtxFqRelSeries}, pr::Tuple{Int, Int} = (5, 2); raw::Bool = false)
  C = G.C
  while precision(C)[1] < pr[1]
    Hecke.MPolyFact.lift_q(C)
  end
  #TODO: truncate precision where necessary, working with an insane
  #      series precision is costly
  while precision(C)[2] < pr[2]
    Hecke.MPolyFact.lift(C)
  end
  rt = [-coeff(x, 0) for x = C.lf[1:C.n]]
  rt = map(y->map_coefficients(x->setprecision(x, pr[1]), setprecision(y, pr[2]), parent = parent(y)), rt)
  if isdefined(G, :rt_num)
    rt = [rt[G.rt_num[i]] for i=1:length(G.rt_num)]
  end
  return rt
end

@doc raw"""
    upper_bound(G::GaloisCtx, f...)

Given a `GaloisCtx` and some multivariate function, upper_bound the image of `f`
upon evaluation at the roots implicit in `G`.

`f` can be
 - a multivariate polynomial or straight-line polynomial (strictly: any object
   allowing `evaluate`
 - `elementary_symmetric` or `power_sum`, in which case more arguments are
   needed: the array with the values and the index.
   `upper_bound(G, power_sum, A, i)` is equivalent to `upper_bound(G, power_sum(A, i))`
   but more efficient.

In every case a univariate polynomial (over the integers) can be added, it
will act as a Tschirnhaus-transformation, ie. the roots (bounds) implicit
in `G` will first be transformed.
"""
function upper_bound end

function upper_bound(G::GaloisCtx, f)
#  @show :eval, f, G.B, degree(G.f)
  return Oscar.evaluate(f, [G.B for i=1:degree(G.f)])
end

function upper_bound(G::GaloisCtx, f, ts::ZZPolyRingElem)
  B = ts(G.B)
  return Oscar.evaluate(f, [B for i=1:degree(G.f)])
end

function upper_bound(G::GaloisCtx, ::typeof(elementary_symmetric), A::Vector, i::Int, ts::ZZPolyRingElem = gen(Oscar.Hecke.Globals.Zx))
  if ts != gen(Hecke.Globals.Zx)
    A = [ts(x) for x = A]
  end
  B = [upper_bound(G, x) for x = A]
  n = length(B)
  b = sort(B)
  return parent(B[1])(binomial(n, i))*prod(b[max(1, n-i+1):end])
end

function upper_bound(G::GaloisCtx, ::typeof(power_sum), A::Vector, i::Int, ts::ZZPolyRingElem = gen(Oscar.Hecke.Globals.Zx))
  if ts != gen(Hecke.Globals.Zx)
    A = [ts(x) for x = A]
  end
  B = [upper_bound(G, x)^i for x = A]
  return sum(B)
end

function upper_bound(G::GaloisCtx, ::typeof(elementary_symmetric), i::Int, ts::ZZPolyRingElem = gen(Oscar.Hecke.Globals.Zx))
  if ts != gen(Hecke.Globals.Zx)
    b = ts(G.B)
  else
    b = G.B
  end
  n = degree(G.f)
  return parent(b)(binomial(n, i))*b^i
end

function upper_bound(G::GaloisCtx, ::typeof(power_sum), i::Int, ts::ZZPolyRingElem = gen(Oscar.Hecke.Globals.Zx))
  if ts != gen(Hecke.Globals.Zx)
    b = ts(G.B)
  else
    b = G.B
  end
  return parent(b)(degree(G.f))*b^i
end

function Hecke.orbit(G::Oscar.PermGroup, f::MPolyRingElem)
  s = Set([f])
  while true
    n = Set(x^g for x = s for g = gens(G))
    sn = length(s)
    union!(s, n)
    if length(s) == sn
      break
    end
  end
  return s
end

function Hecke.evaluate(I::SLPoly, p, a::Vector)
  return evaluate(I, [a[p(i)] for i=1:length(a)])
end

probable_orbit(G::Oscar.PermGroup, f::MPolyRingElem) = orbit(G, f)

"""
`slprogram`s can be compiled into "normal" Julia functions, but there is
some overhead in the compilation itself. By default, apparently nothing is
compiled, so we allow to force this there.

`isPoly` allows the use of inplace operations, as `SLPoly`s result
in programs where intermediate results are used only once.
"""
function compile!(f::SLPoly)
  if !isdefined(f.slprogram, :f)
    Oscar.StraightLinePrograms.compile!(f.slprogram, isPoly = true)
  end
end

#one cannot compare (==) slpoly, no hash either..
#(cannot be done, thus comparison is indirect via evaluation)
#I assume algorithm can be improved (TODO: Max?)
function probable_orbit(G::Oscar.PermGroup, f::SLPoly; limit::Int = typemax(Int))
  n = ngens(parent(f))
  F = GF(next_prime(2^50))
  p = [rand(F) for i=1:n]
  s = [f]
  v = Set([evaluate(f, p)])
  while true
    nw = []
    for g = gens(G)
      for h = s
        z = evaluate(h^g, p)
        if !(z in v)
          push!(nw, h^g)
          push!(v, z)
          if length(s) + length(nw) > limit
            append!(s, nw)
            return s
          end
        end
      end
    end
    if length(nw) == 0
      return s
    end
    append!(s, nw)
  end
end

#TODO:
#- Bessere Abstraktion um mehr Grundkoerper/ Ringe zu erlauben
#- Bessere Teilkpoerper: ich brauche "nur" maximale
#- sanity-checks
#- "datenbank" fuer Beispiele

#a gimmick, not used in galois groups
@doc raw"""
    to_elementary_symmetric(f)

For a multivariate symmetric polynomial `f`, (i.e. `f` is invariant under
permutation of the variables), compute a new polynomial `g` s.th.
`g` evaluated at the elementary symmetric polynomials recovers `f`.

This is using a rather elementary algorithm.

# Examples
We recover the Newton-Girard formulas:
```jldoctest
julia> R, x = polynomial_ring(QQ, 3);

julia> d = power_sum(x, 3)
x1^3 + x2^3 + x3^3

julia> g = to_elementary_symmetric(d)
x1^3 - 3*x1*x2 + 3*x3

julia> evaluate(g, [elementary_symmetric(x, i) for i=1:3])
x1^3 + x2^3 + x3^3
```
"""
function to_elementary_symmetric(f)
  S = parent(f)
  n = ngens(S)
  if n == 1 || is_constant(f)
    return f
  end
  T = polynomial_ring(base_ring(S), n-1)[1]
  g1 = to_elementary_symmetric(evaluate(f, vcat(gens(T), [T(0)])))
  es = [elementary_symmetric(gens(S), i) for i=1:n-1]
  f = f - evaluate(g1, es)
  h = divexact(f, elementary_symmetric(gens(S), n))
  g2 = to_elementary_symmetric(h)
  g1 = evaluate(g1, gens(S)[1:n-1])
  return g1 + gen(S, n)*g2
end

function ^(f::SLPoly, p::Oscar.PermGroupElem)
  #TODO: replace by making the permutation of the input an internal
  #      operation.

  g = gens(parent(f))
  h = typeof(f)[]
  for i=1:ngens(parent(f))
    push!(h, g[p(i)])
  end
  e = evaluate(f, h)
  if typeof(e) != typeof(f)
    @show "bad case"
    return f
  end
  return e
end

@doc raw"""
    isprobably_invariant(g, p) -> Bool

For a multivariate function, mainly an `SLPoly`, test if this is
likely to be invariant under the permutation `p`. Due to the representation
of `SLPoly`s as trees, it is not possible to test this exactly. Instead
`p` is evaluated at random elements in a large finite field.
"""
function isprobably_invariant(g, p)
  R = parent(g)
  k = GF(next_prime(2^20))
  n = ngens(R)
  lp = [rand(k) for i=1:n]
  return evaluate(g, lp) == evaluate(g^p, lp)
end
#TODO: think about the order of arguments!

function isprobably_invariant(p, G::PermGroup)
  R = parent(p)
  k = GF(next_prime(2^20))
  n = ngens(R)
  lp = [rand(k) for i=1:n]
  gp = evaluate(p, lp)
  return all(x->gp == evaluate(p^x, lp), gens(G))
end

function set_orbit(G::PermGroup, H::PermGroup)
  #from Elsenhans
  #https://math.uni-paderborn.de/fileadmin/mathematik/AG-Computeralgebra/inv_transfer_5_homepage.pdf
  #    http://dblp.uni-trier.de/db/journals/jsc/jsc79.html#Elsenhans17
  # https://doi.org/10.1016/j.jsc.2016.02.005

  l = representative.(low_index_subgroup_classes(H, 2*degree(G)^2))
  S, g = slpoly_ring(ZZ, degree(G), cached = false)

  sort!(l, lt = (a,b) -> isless(order(b), order(a)))
  for U = l
    O = orbits(U)
    for o in O
      #TODO: should use orbits of Set(o)...
      f = sum(g[collect(o)])
      oH = probable_orbit(H, f)
      oG = probable_orbit(G, f, limit = length(oH)+5)
      if length(oH) < length(oG)
        for i = 1:length(o)
          I = sum(x^i for x = oH)
          if isprobably_invariant(I, H) &&
             !isprobably_invariant(I, G)
            @vprint :GaloisInvariant 2 "SetOrbit won\n"
            return true, I
          end
        end
        f = prod(g[collect(o)])
        oH = probable_orbit(H, f)
        I = sum(oH)
        if isprobably_invariant(I, H) &&
           !isprobably_invariant(I, G)
          @vprint :GaloisInvariant 2 "SetOrbit won - final attempt\n"
          return true, I
        end
       end
     end
   end
   return false, g[1]
end

@doc raw"""
    invariant(G::PermGroup, H::PermGroup)

For a permutation group `G` and a maximal subgroup `H`, find
a (multivariate (`SLPoly`)) `f` with `G`-stabilizer `H`, i.e. 
`f` is invariant under all permutations in `H`, but not invariant under
any other element of `G`.
"""
function invariant(G::PermGroup, H::PermGroup)
  @vprint :GaloisInvariant 1 "Searching G-relative H-invariant\n"
  @vprint :GaloisInvariant 2 "that is a $G-relative $H-invariant\n"

  S, g = slpoly_ring(ZZ, degree(G), cached = false)

  if is_transitive(G) && !is_transitive(H)
    @vprint :GaloisInvariant 2 "top group transitive, bottom not\n"
    return sum(probable_orbit(H, g[1]))
  end

  if !is_transitive(G) 
    @vprint :GaloisInvariant 2 "both groups are intransitive\n"
    OG = [sort(collect(x)) for x = orbits(G)]
    OH = [sort(collect(x)) for x = orbits(H)]
    d = setdiff(OH, OG)
    if length(d) > 0
      @vprint :GaloisInvariant 2 "groups have different orbits\n"
      return sum(probable_orbit(H, g[d[1][1]]))
    end
    #OH == OG
    for o = OH
      h = action_homomorphism(G, o)
      hG = image(h)[1]
      hH = image(h, H)[1]
      if order(hG) > order(hH)
        @vprint :GaloisInvariant 2 "differ on action on $o, recursing\n"
        @hassert :GaloisInvariant 0 is_maximal_subgroup(hH, hG)
        I = invariant(hG, hH)
        return evaluate(I, g[collect(o)])
      end
    end
    @vprint :GaloisInvariant 2 "going transitive...\n"
    #creating transitive version.
    gs = Oscar.gset(G, [[x[1] for x = OG]])
    os = Oscar.orbits(gs)
    @assert length(os) == 1
    os = os[1]
    h = Oscar.action_homomorphism(os)
    GG = h(G)[1]
    HH = h(H)[1]
    I = invariant(GG, HH)
    ex = 1
    while true
      J = evaluate(I, [sum(g[o])^ex for o = collect(os)])
      if !isprobably_invariant(J, G)
        I = J
        break
      end
      ex += 1
    end
    @hassert :GaloisInvariant 2 isprobably_invariant(I, H)
    @hassert :GaloisInvariant 2 !isprobably_invariant(I, G)

    return I
  end

  if is_primitive(G) && is_primitive(H)
    if isodd(G) && iseven(H)
      @vprint :GaloisInvariant 3 "using sqrt_disc\n"
      return sqrt_disc(g)
    end
    @vtime :GaloisInvariant 2 fl, I = set_orbit(G, H)
    fl && return I
  end

  bG = all_blocks(G)
  bH = all_blocks(H)

  d = setdiff(bH, bG)
  @vprint :GaloisInvariant 2 "Have block systems, they differ at $d\n"
  if length(d) > 0
    @vprint :GaloisInvariant 3  "using F-invar\n"
    return prod(sum(g[b]) for b = block_system(H, d[1]))
  else
    for BB = bG
      #TODO: Max(?): some of this possibly cannot happen anymore
      #      in starting group the sieving has been improved
      #      classical D_sm or similar may no-longer be possible
      B = block_system(H, BB)
      m = length(B)
      l = length(B[1])
      y = [sum(g[b]) for b = B]
      d = [sqrt_disc(g[b]) for b = B]
      D = sqrt_disc(y)

      I = D
      if all(p->isprobably_invariant(I, p), H) &&
         any(p->!isprobably_invariant(I, p), G) #TODO: this can be decided theoretically
        @vprint :GaloisInvariant 3 "using D-invar for $BB\n"
        return I
      end
      I = elementary_symmetric(d, 1)
      if all(p->isprobably_invariant(I, p), H) &&
         any(p->!isprobably_invariant(I, p), G) #TODO: this can be decided theoretically
        @vprint :GaloisInvariant 3 "using s1-invar for $BB\n"
        return I
      end
      I = elementary_symmetric(d, m)
      if all(p->isprobably_invariant(I, p), H) &&
         any(p->!isprobably_invariant(I, p), G) #TODO: this can decided theoretically
        @vprint :GaloisInvariant 3 "using sm-invar for $BB\n"
        return I
      end
      I = elementary_symmetric(d, 2)
      if all(p->isprobably_invariant(I, p), H) &&
         any(p->!isprobably_invariant(I, p), G) #TODO
        @vprint :GaloisInvariant 3 "using s2-invar for $BB\n"
        return I
      end
      I = D*elementary_symmetric(d, m)
      if all(p->isprobably_invariant(I, p), H) &&
         any(p->!isprobably_invariant(I, p), G)
        @vprint :GaloisInvariant 3 "using D sm-invar for $BB\n"
        return I
      end
    end

    for BB = bG
      h = action_on_blocks(G, BB)
      B = block_system(G, BB)
      hG = image(h)[1]
      hH = image(h, H)[1]
      if length(hH) < length(hG)
        y = [sum(g[b]) for b = B]
        J = invariant(hG, hH)
        E = evaluate(J, y)
        @vprint :GaloisInvariant 3 "using E-invar for $BB (4.1.2)\n"
        return E
      end

      sG = stabilizer(G, BB, on_sets)[1]
      sH = stabilizer(H, BB, on_sets)[1]
      hG = action_homomorphism(sG, BB)
      ssG = image(hG, sG)[1]
      ssH = image(hG, sH)[1]
      if length(ssH) < length(ssG)
        J = invariant(ssG, ssH)
        C = left_transversal(H, sH)
        gg = g[BB]
        J = evaluate(J, gg)
        F = sum(J^t for t = C)
        @hassert :GaloisInvariant 1 isprobably_invariant(F, H)
        @hassert :GaloisInvariant 1 !isprobably_invariant(F, G)
        @vprint :GaloisInvariant 3 "using F-invar for $BB (4.1.4)\n"
        return F
      end
    end
  end

  @vprint :GaloisInvariant 3 "no special invar found, resorting to generic\n"

  m = prod(gen(S, i)^i for i=1:degree(G)-1)
  return sum(m^s for s = H)
end

@doc raw"""
    resolvent(C::GaloisCtx, G::PermGroup, U::PermGroup)

Find a `G`-relative `H`-invariant `I` and form the corresponding resolvent polynomial
``\prod (x-I^t)``
where the product runs over all coset-representatives of `G/U`.
"""
function resolvent(C::GaloisCtx, G::PermGroup, U::PermGroup, extra::Int = 5)
  I = invariant(G, U)
  t = right_transversal(G, U)
  n = length(t)
  rt = roots(C)
  #make square-free (in residue field)
  k, mk = residue_field(parent(rt[1]))
  k_rt = map(mk, rt)
  ts = find_transformation(k_rt, I, t)

  B = 2*n*evaluate(I, map(ts, [C.B for i = 1:ngens(parent(I))]))^n
  rt = roots(C, bound_to_precision(C, B, extra))
  rt = map(ts, rt)
  rt = [evaluate(I^s, rt) for s = t]
  pr = copy(rt)
  ps = [isinteger(C, B, sum(pr))[2]]
  while length(ps) < n
    pr .*=  rt
    push!(ps, isinteger(C, B, sum(pr))[2])
  end
  if isa(ps[1], ZZRingElem)
    return Hecke.power_sums_to_polynomial(map(QQFieldElem, ps))
  else
    return Hecke.power_sums_to_polynomial(ps)
  end
end

"""
    minpoly(C::GaloisCtx, I, extra::Int = 5)

Computes the minimal polynomial of `I` evaluated at the roots
stored in `C`.
"""
function Hecke.minpoly(C::GaloisCtx, I, extra::Int = 5)
  O = collect(probable_orbit(C.G, I))
  n = length(O)
  rt = roots(C)
  k = GF(next_prime(10^6))
  r = [rand(k) for i=1:degree(C.f)]
  while length(Set(r)) < length(r)
    r = [rand(k) for i=r]
  end
  ts = find_transformation(r, O)

  B = 2*n*evaluate(I, map(ts, [C.B for i = 1:ngens(parent(I))]))^n
  rt = roots(C, bound_to_precision(C, B, extra))
  rt = map(ts, rt)
  rt = [evaluate(i, rt) for i = O]
  pr = copy(rt)
  ps = [isinteger(C, B, sum(pr))[2]]
  while length(ps) < n
    pr .*=  rt
    push!(ps, isinteger(C, B, sum(pr))[2])
  end
  if isa(ps[1], ZZRingElem)
    return Hecke.power_sums_to_polynomial(map(QQFieldElem, ps))
  else
    return Hecke.power_sums_to_polynomial(ps)
  end
end

Hecke.minpoly(R::ZZPolyRing, C::GaloisCtx, I, extra::Int = 5) = Hecke.minpoly(C, I, extra = extra)(gen(R))
Hecke.minpoly(R::QQPolyRing, C::GaloisCtx, I, extra::Int = 5) = Hecke.minpoly(C, I, extra = extra)(gen(R))

struct GroupFilter
  f::Vector{Tuple{Function, String}}
  GroupFilter() = new([(x->true, "")])
end

#TODO: Max: rank the filter function by cost, cheapest first...
#      possibly needs a third field in the structure
function (F::GroupFilter)(G::PermGroup)
  do_print = Hecke.get_verbosity_level(:GaloisGroup) >= 1
  do_all = Hecke.get_verbosity_level(:GaloisGroup) >= 2

  res = true
  for (x, s) = F.f
    fl = x(G)
    fl && continue
    if do_all
      if res 
        println("group rejected by:\n\t- $s")
      elseif length(s) > 0
        println("\t- $s")
      end
      res = false
    elseif do_print
      println("group rejected by: $s")
      return false
    else
      return false
    end
  end
  return res
end

function Base.push!(G::GroupFilter, F::Function, s::String="")
  push!(G.f, (F, s))
end

@doc raw"""
# `DescentEnv`

Given a permutation group `G`, known to contain the Galois group of some polynomial,
the key step in Stauduhar's method is to systematically test all maximal subgroups
of `G`. This structure is used to organize this process. On creation, all the
maximal subgroups are computed and possibly filtered (if the polynomial
is irreducible, only transitive groups are needed, if the stem field if primitive,
then so is the group...). The filtering is done via a `GroupFilter`.

Furthermore we may compute (and store)
 - a `G`-relative `H`-invariant for any subgroup (`I`)
 - a list of Tschirnhausen transformations used (`T`)
 - a list (`l`) marking which subgroups we have dealt with (`1` means group has been
 processed completely)
"""
mutable struct DescentEnv
  G::PermGroup
  s::Vector{PermGroup}
  I::Dict{Int, SLPoly}
  T::Dict{Int, Vector{ZZPolyRingElem}}
  l::Vector{Int}
  rng::AbstractRNG
  #the coset reps need to be cached as well
  #a "work limit" on the "invariant" function
  #a more select choice of group....

  function DescentEnv(G::PermGroup, f::GroupFilter = GroupFilter())
    s = map(representative, maximal_subgroup_classes(G))
    r = new()
    r.G = G
    @vprint :GaloisGroup 1 "starting with $(length(s)) maximal subgroup classes\n"
    r.s = filter(f, s)
    @vprint :GaloisGroup 1 "filtering down to $(length(r.s)) maximal subgroup classes\n"
    r.I = Dict{Int, SLPoly}()
    r.T = Dict{Int, Vector{ZZPolyRingElem}}()
    r.l = zeros(Int, length(r.s))
    r.rng = MersenneTwister(1)
    return r
  end
end

#TODO: instead of iterating, use the cheapest invariant
#      use invariant with a worklevel to limit costs in finding them
#      use the V_4 (index 2) stuff when necessary
function Base.iterate(D::DescentEnv, i::Int=1)
  #r.l[i] == 1: done,  0:: nothing
  all(isone, D.l) && return nothing

  while i<=length(D.s) && D.l[i] == 1
    i += 1
  end

  local I, TS
  if haskey(D.I, i)
    I = D.I[i]
  else
    I = D.I[i] = invariant(D.G, D.s[i])
  end

  if haskey(D.T, i)
    T = D.T[i]
    ts = rand(D.rng, Hecke.Globals.Zx, 2:degree(D.G), -4:4)
    while degree(ts) < 2 || ts in T
      ts = rand(D.rng, Hecke.Globals.Zx, 2:degree(D.G), -4:4)
    end
    push!(D.T[i], ts)
    TS = ts
  else
    D.T[i] = [gen(Hecke.Globals.Zx)]
    TS = D.T[i][1]
  end

  return (D.s[i], I, TS), i
end

function Base.push!(D::DescentEnv, i::Int)
  D.l[i] = 0
end

function Base.pop!(D::DescentEnv, i::Int)
  D.l[i] = 1
end

function sum_orbits(K, Qt_to_G, r)
  @vprint :GaloisGroup 1 "group is primitive (no subfields), trying operations on pairs\n"

  #starting group as the stabilizer of the factorisation of the 2-sum poly,
  #the polynomial with roots r_i + r_j, i<j
  #we need this square-free, so we compute this over the finite field
  #actually, we only check that the roots over the finite field are distinct.
  #if not: transform (using ts) and try again.
  
  k = parent(r[1])
  m = Dict{typeof(r[1]), Vector{Int}}()  #the vector has length 2 always
  local ts = gen(Hecke.Globals.Zx)
  c = r
  while true
    for i=1:length(r)-1
      for j=i+1:length(r)
        m[c[i]+c[j]] = [i,j]
      end
    end
    if length(keys(m)) < binomial(degree(K), 2)
      @vprint :GaloisGroup 2 " 2-sum: found duplicate, transforming...\n"
      while true
        ts = rand(Hecke.Globals.Zx, 2:degree(K), -4:4)
        if degree(ts) > 1
          break
        end
      end
      c = map(ts, r)
      empty!(m)
    else
      break
    end
  end

  @vprint :GaloisGroup 1 "have everything, now getting the 2-sum poly\n"
  if gen(Hecke.Globals.Zx) == ts
    @vtime :GaloisGroup 2 g = msum_poly(defining_polynomial(K), 2) #if f has roots a_i, then g has roots a_i+a_j, if G is 2-transitive, this is pointless.
  else
    @vtime :GaloisGroup 2 g = msum_poly(minpoly(ts(gen(K))), 2)
  end
  @vprint :GaloisGroup 1 "... factoring...\n"
  @vtime :GaloisGroup 2 fg  = factor(g)
  @assert all(isone, values(fg.fac))

  O = []
  if isa(r[1], AcbFieldElem)
    mm = collect(m)
  end
  for f = keys(fg.fac)
    r = roots(map_coefficients(Qt_to_G, f))
    if isa(r[1], AcbFieldElem)
      push!(O, [mm[argmin(map(x->abs(x[1]-y), mm))][2] for y = r])
    else
      push!(O, [m[x] for x = r])
    end
  end
  @vprint :GaloisGroup 2 "partitions: $O\n"
  return O
end

"""
Given 2 normal subgroups of index 2, there is a 3 index 2 subgroup
automatically, given that the quotient ov G by the intersection
is a V_4 group.

This computes (in a dumb way) the 3rd group
"""
function index2sum(G::PermGroup, H::PermGroup, V::PermGroup)
  @assert index(G, H) == index(G, V) == 2
  U = intersect(H, V)[1]
  @assert index(G, U) == 4
  t = right_transversal(G, U)
  s = [x for x = t if !(x in H || x in V)]
  @assert length(s) == 1
  return sub(G, vcat(gens(U), s))[1]
end

"""
Maps elements of the base ring into the completion (the splitting field).
FOr computations over QQ this does not matter as QQ just coerces into
any q/p-adic field.

FOr computations over number fields this is not true, here an explciti map
is needed.

TODO: check and add a precision parameter
TODO: implement also for Q(t) (even though not used)
"""
function map_coeff(C::GaloisCtx{Hecke.qAdicRootCtx}, x)
  return (C.C.Q[1])(x)
end

@doc raw"""
    starting_group(GC::GaloisCtx, K::SimpleNumberField) 

Finds a _starting group_, that is a group `G` as a subgroup of the
symmetric group acting on the roots in the explicit ordering in the
galois context.

If the field is imprimitive, the group is derived from the subfields, otherwise
the factorisation of the 2-sum polynomial is used.

Returns a triple: 
 - the group
 - a filter for all groups that can occur
 - a permutation representing the operation of the Frobenius automorphism
"""
function starting_group(GC::GaloisCtx, K::T; useSubfields::Bool = true) where T <: SimpleNumField
  c = roots(GC, 5, raw = true)

  @vprint :GaloisGroup 1 "computing starting group (upper bound for Galois group)\n"

  if useSubfields
    @vprint :GaloisGroup 1 "computing principal subfields ...\n"
    @vtime :GaloisGroup 2 S = Hecke.principal_subfields(K)
  else
    S = []
  end

  #compute the block system for all subfields...
  bs = Vector{Vector{Int}}[]
  fld = Dict{Vector{Vector{Int}}, Any}()
  for (s, ms) = S
    if degree(s) == degree(K) || degree(s) == 1
      continue
    end
    pr = 5
    local v
    #TODO: maybe also use cycle_types for lower bounds?
    while true
      c = roots(GC, pr, raw = true)
      g = ms(gen(s))
      gg = map_coefficients(x->map_coeff(GC, x), parent(K.pol)(g))
      d = map(gg, c)
      v = Hecke.MPolyFact.block_system(d)
      if any(x->length(x) != div(degree(K), degree(s)), v)
        pr *= 2
        if pr > 100
          error("too bad")
        end
      else
        break
      end
    end
    #each subfield defines a block system, each block system results in
    #a wreath product representing the largest subgroup of Sn with this
    #block system: Sl wr Sk
    #However, this wreath prodct has automatically 3 maximal subgroups
    #of index 2 (Eichenlaub Prop 3, p 28 or Geissler Satz 6.8).
    # - Sl wr Ak: possible iff disc(subfield) is a square, so
    #   we can test this and then reduce the starting group 
    #      or
    #   add this as an exclusion: we don't need to test subgroups of this
    #   any more
    # - one of them should be (Sl wr Sk) meet A(lk))
    #   disc(K) is a square
    # - see below in starting_group
    sort!(v, lt = (x,y) -> minimum(x) < minimum(y))
    push!(bs, v)
    fld[v] = (s, ms)
  end

  F = GroupFilter()
  push!(F, is_transitive, "not transitive") #poly is defining number field, hence irreducible


  #selecting maximal block systems only...
  if length(bs) > 1
    @vprint :GaloisGroup 2 "group will have (all) block systems: $([x[1] for x = bs])\n"
    L = POSet(bs, (x,y) -> issubset(x[1], y[1]) || issubset(y[1], x[1]),
                  (x,y) -> issubset(x[1], y[1]) - issubset(y[1], x[1]))
    bs = map(Oscar.data, minimal_elements(L))
    @vprint :GaloisGroup 1 "group will have (maximal) block systems: $([x[1] for x = bs])\n"
  end
  GC.start = (1, bs)

  G = symmetric_group(degree(K))
  if degree(K) == 1
    return G, F, one(G)
  end
  S = G

  #we only need groups that allow the subfield to keep the
  #"correct" galois group. Up to isomorphism, the subfield galois group
  #is the operation on the block-system (permuting of blocks)
  for v = bs
    s, ms = fld[v] #it appears that principal_subfields are maximal..
                   #(all others are via field intersection, hence smaller)
    if degree(s) > 2
      #no point in checking C_2
      local gc, bs
      # TODO: infer the subfields of s from the ones in K
      #       the poset should make this trivial.
      @vprint :GaloisGroup 1 "computing group of subfield of degree $(degree(s))\n"
      local g, gc
      if isa(parent(c[1]), FlintLocalField) && 
         valuation(discriminant(s), GC.prime) > 0
        throw(Hecke.BadPrime(GC.prime))
      end
      g, gc = try
          Hecke.pushindent()
          if isa(c[1], NumFieldElem)
            galois_group(s, algorithm = :Symbolic, field = parent(c[1]))
          elseif isa(c[1], FlintLocalFieldElem)
            galois_group(s, prime = GC.prime)
          else
            galois_group(s, algorithm = :Complex)
          end
        catch e  #can "legally" be BadPrime, will be dealt with one level up
          Hecke.popindent()
          rethrow()
        end
      Hecke.popindent()

      i = order(g)

      pr = 5
      r = roots(gc, pr, raw = true)
      R = roots(GC, pr, raw = true)

      gg = map_coefficients(x->map_coeff(GC, x), parent(K.pol)(ms(gen(s))))
      d = map(gg, R)
      if isa(r[1], AbsSimpleNumFieldElem)
        @assert parent(r[1]) == parent(R[1])
        f = _F = parent(r[1])
        mf = mF = x->x
        mfF = x->x
      elseif isa(r[1], AcbFieldElem)
        f = _F = parent(r[1])
        mf = mF = x->x
        mfF = x->x
      else
        f, mf = residue_field(parent(r[1]))
        _F, mF = residue_field(parent(R[1]))
        mfF = Hecke.find_morphism(f, _F)
      end
      #we should have
      # - d == r (in the appropriate setting)
      # - order(Set(map(mF, d))) == order(Set(map(mf, r))) == degree(s)
      # - Set(map(mfF, map(mf, r))) == Set(map(mF, d))
      # - d is partitioned
      # and g (the small Galois group) operates on the blocks in d
      # Now we need to create the wreath_product(sym(deg(K)), g)
      # that produces THIS block system
      # The ambiguity of the different possible embeddings (mfF)
      # should be harmless as they are linked by the Frobenius
      # and Frobenius in g, so g^Frobenius == g....
      r = map(mf, r)
      r = map(mfF, r)
      @assert parent(R[1]) == parent(d[1])
      d = map(mF, d)
      R = map(mF, R)
      if length(Set(r)) == length(r) &&
         length(Set(d)) == divexact(length(d), length(r))
         length(Set(R)) == length(R) #no problems with precision, 
                                     #we can proceed.
        @assert Set(r) == Set(d)
        #assuming wreath_product(A, B) has block system [[1..l], [l+1..2l]...]
        # the re-ordering is easy: W^s for s = vcat(bs)
        bs = map(x->findall(isequal(x), d), r)
        @assert all(x->length(x) == length(bs[1]), bs)
        W = PermGroup(wreath_product(symmetric_group(length(bs[2])), g))
        #should have the block system as above..
        W = W^symmetric_group(degree(W))(vcat(bs...))
        G = intersect(G, W)[1]
      else
        @vprint :GaloisGroup 2 "precision problem in wreath product, using size only\n"
      end
      let _v = v, _i = i 
        push!(F, x->order(image(action_on_blocks(x, _v[1]))[1]) % _i == 0 , "subfield for $(_v[1]) would have wrong group (too small)")
      end
    end
  end

  #TODO: make generic!!!
  if isa(c[1], AcbFieldElem)
    d = map(conj, c)
    si = [argmin(map(y->abs(y-x), c)) for x = d]
  elseif isa(c[1], AbsSimpleNumFieldElem) #.. and use automorphism
    si = collect(1:length(c))
  else
    d = map(frobenius, c)
    si = [findfirst(y->y==x, c) for x = d]
  end

  @vprint :GaloisGroup 1 "found Frobenius element: $si\n"

  in_br = false
  if is_square(discriminant(K))
    G = intersect(G, alternating_group(degree(K)))[1]
    in_br = true
  else
    push!(F, isodd, "has to be odd")
  end

  for b = bs
    #each block system b corresponds to a subfield and induces an upper
    #limit of i
    #  G =  Sym(k) wr Sym(l) 
    #on the Galois group
    #This wreath product has always (at least) 3 index 2 subgroups:
    #  G meet Alt(k*l) 
    #  Sym(k) wr Alt(l)
    #  index2sum of the above
    #We are a subgroup of the 1st iff disc(K) is a square (in_br = true)
    #                         2nd iff disc(k) is a square (in_ar = true)
    #                         3rd iff disc(k)*disc(K) is a square
    #This is used in both directions
    # - in an intersection to limit the starting group (upper bound)
    #   (true case)
    # - in the filter to avoid unnecessary tests
    #   (false case)
    # disc(k) is up to an unknown Tschirnhaus transformation exactly the
    # D(y_1, .. y_k) invar in Geisser/Eichenlaub (case 2)
    # and the product then is one for the 3rd case.
    s = S(vcat(b...))

    wr = PermGroup(wreath_product(symmetric_group(length(b[1])), symmetric_group(length(b))))
    br = intersect(wr, alternating_group(degree(K)))[1]
    wr = wr^s
    br = br^s

    G = intersect(G, wr)[1]

    can_use_wr = length(b) > 2
    #for length(b) == 2, the bottom field is quadratic and thus always S2
    #the wreath product would not be transitive here
    #TODO: the new stuff above should guarantee that the induced action
    #  on a block system can never be too large:
    #    (sym(l) wr gal(k) is constructed above
    #    does this part still make sense????
    if can_use_wr
      ar = PermGroup(wreath_product(symmetric_group(length(b[1])), alternating_group(length(b))))
      ar = ar^s
      cr = index2sum(wr, ar, br)
    end


    in_ar = false
    if is_square(discriminant(fld[b][1]))
      in_ar = true
    end
    #in_ar == true iff galois <= ar
    #in_br == true iff galois <= br
    in_cr = is_square(discriminant(K)*discriminant(fld[b][1]))

    if can_use_wr
      if in_ar 
        G = intersect(G, ar)[1]
      else
        let ar = ar 
          push!(F, x->!is_subset(x, ar), "subfield is even/odd")
        end
      end
    end
    if in_br
      #G = intersect(G, br)[1] #already done above
    else
      #push!(F, x->!is_subset(x, br)) #already done above
    end
    if can_use_wr
      if in_cr
        G = intersect(G, cr)[1]
      else
        let cr = cr
          push!(F, x->!is_subset(x, cr), "third subgroup of wreath product")
        end
      end
    end
  end

  if length(bs) == 0 #primitive case: no subfields, no blocks, primitive group!
    push!(F, is_primitive, "primitivity")
    pc = parent(c[1])
    if isa(pc, NumField)
      k = pc
      mk = x->x
    elseif isa(c[1], AcbFieldElem)
      k = pc
      mk = x->x
    else
      k, mk = residue_field(pc)
    end
    O = sum_orbits(K, x->mk(pc(map_coeff(GC, x))), map(mk, c))
    GC.start = (2, O)
    
    #the factors define a partitioning of pairs, the stabilizer of this
    #partition is the largest possible group...
    #code from Max...

    #TODO: wrap this properly
    if hasproperty(GAP.Globals, :ConStabilize)
      tmp = [GAP.Globals.ConStabilize(GAP.Obj(sort(o), recursive=true), GAP.Globals.OnSetsSets) for o in O]
      H = GAP.Globals.Solve(GAP.Obj(vcat(GAP.Globals.ConInGroup(G.X), tmp)))
      G = Oscar._as_subgroup(G, H)[1]
    else
      #TODO: fallback if ferret wasn't loaded for some reason; this should be removed
      # once we have GAP_pkg_ferret
      G = intersect([stabilizer(G, sort(o), on_sets_sets)[1] for o=O]...)[1]
    end
  end

  si = G(si)

  return G, F, si
end

@doc raw"""
    find_prime(f::QQPolyRingElem, extra::Int = 5; prime::Int = 0, pStart::Int = 2*degree(f)) -> p, c

Tries to find a useful prime for the computation of Galois group. Useful means
  - degree of the splitting field not too small
  - degree of the splitting field not too large
Starts searching at `pStart`, returns the prime as well as the cycle types identified.

If `prime` is given, no search is performed.
"""
function find_prime(f::QQPolyRingElem, extra::Int = 5; prime::Int = 0, pStart::Int = 2*degree(f), filter_prime = x->true, filter_pattern = x->true)
  if prime != 0
    p = prime
    lf = factor(GF(p), f)
    return p, Set([CycleType(map(degree, collect(keys(lf.fac))))])
  end
  if pStart < 0
    error("should no longer happen")
    p = -pStart
    lf = factor(GF(p), f)
    return p, Set([CycleType(map(degree, collect(keys(lf.fac))))])
  end

  d_max = degree(f)^2
  d_min = min(2, d_max)
  p_best = 1
  cnt = 5
  ct = Set{Vector{Int}}()
  ps = Vector{Tuple{Int, Int}}()

  #find a q-adic splitting field of "good degree":
  # - too small, then the Frobenius automorphisms is not containing lots of
  #   information
  # - too large, all is slow.
  # careful: group could be (C_2)^n hence d_min might be small...
  no_p = 0
  for p = Hecke.PrimesSet(pStart, -1)
    filter_prime(p) || continue
    k = GF(p)
    if k(leading_coefficient(f)) == 0
      continue
    end
    lf = factor(GF(p), f)
    if any(x->x>1, values(lf.fac))
      continue
    end
    filter_pattern(lf) || continue
    no_p +=1 
    push!(ct, sort(map(degree, collect(keys(lf.fac))))) # ct = cycle types as vector of cycle lengths
    d = lcm([degree(x) for x = keys(lf.fac)])
    if d <= d_max 
      if d >= d_min
        push!(ps, (p, d))
      else
        no_p += 1
        if no_p > degree(f)
          d_min = 1
        end
      end
    end
    if length(ps) > 2*degree(f)
      break
    end
  end

  @vprint :GaloisGroup 2 "possible primes and degrees: $ps\n"

  #want: degree > degree K/2 - which might not be possible.
  l = [x for x = ps if 2*x[2] <= degree(f)]
  if length(l) > 0
    pd = l[argmax([x[2] for x = l])]
  else
    pd = ps[argmin([x[2] for x = ps])]
  end

  p = p_best = pd[1]
  d_max = pd[2]

  @vprint :GaloisGroup 1 "using prime $p_best with degree $d_max\n"
  @vprint :GaloisGroup 2 "and cycle types $ct\n"
  return p_best, Set(CycleType(x) for x = ct)
end

#given cycle types, try to find a divisor (large) of the transitive
#group. Used for sieving.
# Elsenhans-Klueners:
# Computing subfields of number fields and
# applications to Galois group computations
function order_from_shape(ct::Set{CycleType}, n)
  o1 = ZZRingElem(1)
  #Lemma 15:
  o  = n * reduce(lcm, ZZRingElem[divexact(
           reduce(lcm, ZZRingElem[l for (l,v) = x]),
           reduce(gcd, ZZRingElem[l for (l,v) = x])) for x = ct if length(x) > 0], 
           init = ZZRingElem(1))
  #Lemma 16:         
  for p = PrimesSet(1, n)
    #find cycletypes that would correpond to elt of order a multiple of p
    #the exact order is lcm(x) for x in ct
    ct_p = [x for x = ct if any(t->t[1] % p == 0, x)]
    #so ct_t contains all cycle types belonging to elements where the
    #order is divisible by p
    length(ct_p) == 0 && continue
    #now x in ct_p is is the cycle type of an element where the order is a multiple
    #of p. I need the type of x^f (which is of maximal order p-power order)
    #e.g. a cycle of length 6 will produce 3 of length 2
    ct_pp = Set{CycleType}()
    for x = ct_p
      y = x^Int(ppio(order(x), ZZRingElem(p))[2])
      push!(ct_pp, y^y.s[1][1])
    end
    # Now everything has a fixed point, so we can bound the
    # size of Stab_G(1) from below (we're transitive, so the fixed point
    # is 1)
    # the plan is to take to 2 different largest ones

    mu = [(y.s[end][1], y.s[end][1]*y.s[end][2]) for y = ct_pp]
    sort!(mu, lt = (a,b) -> a[1] < b[1])
    if length(mu) == 1
      p_part = mu[1][1]
      o1 *= p_part
    else
      p_part = maximum([x[1]*y[1] for x = mu for y = mu if x[2] != y[2]], init = maximum(x[1] for x = mu))
      o1 *= p_part
    end
  end
  return lcm(n*o1, o)
end


@doc raw"""
    primitive_by_shape(ct::Set{CycleType}, n::Int)

Return `true` if a transitive group $G \leq Sym(n)$ containing permutations
matching the cycle types in `ct` must be primitive, otherwise return `false`.
"""
function primitive_by_shape(ct::Set{CycleType}, n::Int)

  # degree <= 3: always A_n or S_n
  n <= 3 && return true

  # if n is a prime then any transitive action is primitive
  is_prime(n) && return true

  # Now use the following observation:
  #
  # Lemma: Let G \leq Sym(n) be transitive. If there is \pi \in G with a
  # non-trivial cycle whose length p is prime and satisfies n/2 < p < n,
  # then G is primitive.
  #
  # Proof: If B is a block of length b, and n = bk, then G is a subgroup of
  # S_b \wr S_k, which has size (b!)^k * k!. Suppose 1<b<n holds, then this
  # order is not divisible by p, as p > n/2. But the order of \pi is a
  # multiple of p. Contradiction!
  for p in PrimesSet(div(n, 2)+1, n-1)
    if any(y->any(x->x[1] == p, y), ct)
      return true
    end
  end

  # degree 9: also primitive if it contains a double-4-cycle
  # (as the only possible non-trivial block size is 3, and it is easy
  # to see that this is impossible when there is a double-4-cycle.
  n == 9 && return any(x -> (4 => 2) in x, ct)

  # degree 10: also primitive if it contains a triple-3-cycle
  # (as the only possible non-trivial block size are 2 and 5, and it is easy
  # to see that this is impossible when there is a triple-3-cycle.
  n == 10 && return any(x -> (3 => 3) in x, ct)

  # degree 15: also primitive if it contains a double-7-cycle
  # (but this still doesn't catch all cases)
  n == 15 && return any(x -> (7 => 2) in x, ct)

  # TODO: if desired, we could add more degrees here if we wanted, possibly
  # even in a semi-generic way. One could also write a program which uses the
  # databases of transitive and primitive groups to identify additional
  # sufficient conditions.

  return false
end


@doc raw"""
    an_sn_by_shape(ct::Set{CycleType}, n::Int)

Return `true` if a transitive group $G \leq Sym(n)$ containing permutations
matching the cycle types in `ct` must contain $Alt(n)$, otherwise return `false`.
"""
function an_sn_by_shape(ct::Set{CycleType}, n::Int)
  # handle small degrees first

  # degree <= 3: always A_n or S_n
  n <= 3 && return true

  # degree 4 and 5: must contain a 3 cycle
  (n == 4 || n == 5) && return any(x -> (3 => 1) in x, ct)

  # degree 6: must contain a 3 cycle and a 5 cycle
  n == 6 && return any(x -> (3 => 1) in x, ct) && any(x -> (5 => 1) in x, ct)

  # degree 7: must contain a 5 cycle OR contain "many" cycle types
  n == 7 && return length(ct) > 5 || any(x -> (5 => 1) in x, ct)

  # In general, we use a result due to Jordan: Let G <= S_n be transitive, and let
  # p be a prime satisfying n/2 < p < n − 2. If G contains an element that has
  # a cycle of length p then G is a giant (i.e., is A_n or S_n).
  #
  # See e.g. Corollary 10.2.2 on page 225 in [Ser03].
  #
  # Note that this condition is sufficient but not necessary. Indeed for n <= 7
  # no such p exists, which is why we have to treat these degrees separately above.
  lo = div(n, 2)+1
  hi = n-3
  return any(x->any(y -> lo <= y[1] <= hi && isodd(y[1]) && is_prime(y[1]), x.s), ct)
end

#
#TODO
# - subfield lattice directly
# - subfield to block system outsource (also for lattice)
# - more invars
# - re-arrange to make cheap things first, expensive later
# - proof, live and later
# - tschirnhausen transformation: slowly increasing degree and coeffs.
# - for larger degrees: improve starting group by doing Galois groups of subfields
# - more base rings
# - applications: subfields of splitting field (done), towers, solvability by radicals
@doc raw"""
    galois_group(K::AbsSimpleNumField, extra::Int = 5; useSubfields::Bool = true, pStart::Int = 2*degree(K)) -> PermGroup, GaloisCtx

Computes the Galois group of the splitting field of the defining polynomial of `K`.
Currently the polynomial needs to be monic.

The group is returned as an explicit permutation group permuting the roots as contained
in the context object (the 2nd return value). The roots live in a suitable unramifed
extension of the p-adics.

# Examples
```jldoctest
julia> K, a = cyclotomic_field(5);

julia> G, C = galois_group(K)
(Permutation group of degree 4 and order 4, Galois context for x^4 + x^3 + x^2 + x + 1 and prime 19)

julia> describe(G)
"C4"

julia> roots(C, 2)
4-element Vector{QadicFieldElem}:
 (4*19^0 + 2*19^1 + O(19^2))*a + 5*19^0 + 9*19^1 + O(19^2)
 (15*19^0 + 16*19^1 + O(19^2))*a + 9*19^0 + 7*19^1 + O(19^2)
 (18*19^0 + 18*19^1 + O(19^2))*a + 12*19^0 + O(19^2)
 (19^0 + O(19^2))*a + 11*19^0 + 19^1 + O(19^2)
```
"""
function galois_group(K::AbsSimpleNumField, extra::Int = 5; 
  useSubfields::Bool = true, 
  pStart::Int = 2*degree(K), 
  prime::Int = 0, 
  do_shape::Bool = true,
  algorithm::Symbol=:pAdic, 
  field::Union{Nothing, AbsSimpleNumField} = nothing)

  @assert algorithm in [:pAdic, :Complex, :Symbolic]

  if do_shape
    p, ct = find_prime(K.pol, pStart)
  else
    ct = Set{CycleType}()
    @assert prime != 0
  end

  if prime != 0
    p = prime
  end

  # TODO: otherwise, try to detect here if we are primitive or even 2-transitive
  #       primitive: done.
  #       2-transitive: ????
  # TODO: sieve (filter) by 
  #       - cycle_structure (done)
  #       - primitive/ transitive/ 2-transitive/ ... (partly done)
  #       - blocks (not sure: we're using the subfield group and wreath prods)

  # TODO: this too needs to be made either generic or replicated in the 
  #       relative case
  #       or the RootCtx needs to learn to deal with bad primes
 
  while true
    if algorithm == :pAdic
      GC = GaloisCtx(Hecke.Globals.Zx(numerator(K.pol)), p)
    elseif algorithm == :Complex
      GC = GaloisCtx(Hecke.Globals.Zx(numerator(K.pol)), AcbField(20))
    elseif algorithm == :Symbolic
      GC = GaloisCtx(Hecke.Globals.Zx(numerator(K.pol)), field)
    else
      error("wrong algorithm used")
    end
    r = roots(GC, 5)
    if length(r) < degree(K) 
      @vprint :GaloisGroup 1 "in recursion: bad prime detected\n"
      algorithm != :pAdic || error("internal error aborting")
      prime == 0 || throw(Hecke.BadPrime(p))
      p, _ = find_prime(K.pol, pStart = p+1)
      continue
    end

    if an_sn_by_shape(ct, degree(K))
      @vprint :GaloisGroup 1 "An/Sn by cycle type\n"
      if is_square(discriminant(K))
        G = alternating_group(degree(K))
      else
        G = symmetric_group(degree(K))
      end
      GC.G = G
      return G, GC
    end

    G, F, si = try
      starting_group(GC, K, useSubfields = useSubfields && !primitive_by_shape(ct, degree(K)))
    catch e
      if isa(e, Hecke.BadPrime)
        @vprint :GaloisGroup 1 "in recursion: bad prime detected\n"
        prime == 0 || rethrow()
        p, _ = find_prime(K.pol, pStart = p+1)
        continue
      else
        rethrow(e)
      end
    end
    # the filter will filter at least on 
    #   size (have a lower bound from cycles)
    #   even/odd
    #   some subgroups of the wreath products
    if degree(K) == 1
      GC.G = G
      return G, GC
    end

    os = order_from_shape(ct, degree(K))
    @vprint :GaloisGroup 2 "using lower bound $os for order to sieve\n"
    push!(F, x->order(x) % os == 0, "group too small")
    if degree(K) < 20 #arbitrary....
      push!(F, x->issubset(ct, cycle_structures(x)), "by cycle structures")
    end

    # TODO: here we know if we are primitive; can we detect 2-transitive (inside starting_group)?

    return descent(GC, G, F, si, extra = extra)
  end
end

@doc raw"""
    descent(GC::GaloisCtx, G::PermGroup, F::GroupFilter, si::PermGroupElem; grp_id = transitive_group_identification, extra::Int = 5)

Performs a generic Stauduhar descent: starting with the group `G` that needs to be a 
supergroup of the Galois group, operating on the roots in `GC`, the context object.

The groups are filtered by `F` and the result needs to contain the permutation `si`.
For verbose output, the groups are printed through `grp_id`.
"""
function descent(GC::GaloisCtx, G::PermGroup, F::GroupFilter, si::PermGroupElem; grp_id = transitive_group_identification, extra::Int = 5)
  @vprint :GaloisGroup 2 "Have starting group with id $(grp_id(G))\n"

  n = degree(GC.f)

  nG = order(G)
  while true
    D = DescentEnv(G, F)
    @vprint :GaloisGroup 2 "found $(length(D.s)) maximal subgroups\n"

    while true
      d = iterate(D)
      if d === nothing
        break
      end
      s, I, ts = d[1]

      @vprint :GaloisGroup 2 "testing descent $(grp_id(G)) -> $(grp_id(s)) of index $(index(G, s))\n"

      @hassert :GaloisInvariant 1 all(x->isprobably_invariant(I, x), s)
      @hassert :GaloisInvariant 1 any(x->!isprobably_invariant(I, x), G)

      B = upper_bound(GC, I, ts)
      @vprint :GaloisGroup 2 "invariant uses $(cost(I, ts)) multiplications\n"
      @vprint :GaloisGroup 2 "of maximal degree $(total_degree(I)*degree(ts))\n"
      @vprint :GaloisGroup 2 "root upper_bound: $(value(B))\n"

      c = roots(GC, bound_to_precision(GC, B, extra))
      c = map(ts, c)

      local fd
      cnt = 0

      cs = Set{typeof(c[1])}()
      fd = []

      local lt
      if index(G, s) < 100
        @vtime :GaloisGroup 2 lt = right_transversal(G, s)
      elseif is_normalized_by(s, G)
        lt = [one(G)] # I don't know how to get the identity
      else
        @vtime :GaloisGroup 2 lt = short_right_transversal(G, s, si)
        @hassert :GaloisInvariant 2 all(x->si in s^x, lt)
      end
      #x^7-2 triggers wrong descent, occaisonly, if lt is only 10 elems long
      while length(lt) < 20 && length(lt) < index(G, s)
        r = rand(G)
        if all(x->right_coset(s, r) != right_coset(s, x), lt)
          push!(lt, r)
        end
      end

      compile!(I)
      for t = lt
        e = evaluate(I, t, c)
        if typeof(e) == AcbFieldElem && any(x->abs(e-x) < 1e-10, cs)
          @vprint :GaloisGroup 2 " evaluation found duplicate, transforming...\n"
          push!(D, d[2])
          break
        end
        if e in cs
          @vprint :GaloisGroup 2 " evaluation found duplicate, transforming...\n"
          push!(D, d[2])
          break
        end

        push!(cs, e)
        fl, l = isinteger(GC, B, e)
        if fl
          @vprint :GaloisGroup 2 "found descent at $t and value $l\n"
          push!(fd, t)
        end
      end
      if length(cs) == length(lt)
        pop!(D, d[2])
        if length(fd) > 0
          @vprint :GaloisGroup 2 "descending via $fd\n"
        else
          @vprint :GaloisGroup 2 "no descent\n"
        end
        if length(fd)>0
          push!(GC.chn, (G, I, ts, fd))
          G = intersect([s^x for x = fd]...)[1]
          @assert length(fd) > 1 || order(G) == order(s)
          @assert order(G) < nG
          @vprint :GaloisGroup 2 "descending to $(grp_id(G))\n"
          
          if order(G) == n
            GC.G = G
            return G, GC
          end
          break
        end
      end
    end
    if nG == order(G)
      GC.G = G
      return G, GC
    end
    nG = order(G)
  end      
  GC.G = G
  return G, GC
end

"""
    isinteger(C::GaloisCtx, B::BoundRingElem, v)

For an element `v` representing an integral polynomial evaluated at the
roots stored in `C`, known to be bounded from above by `B`, either return
`true` and an explicit (algebraic) integer in the base ring of the context or
return `false`.
"""
function isinteger end

function isinteger(GC::GaloisCtx{Hecke.qAdicRootCtx}, B::BoundRingElem{ZZRingElem}, e)
  p = GC.C.p
  if e.length<2
    l = coeff(e, 0)
    lz = lift(l)
    lz = Hecke.mod_sym(lz, ZZRingElem(p)^precision(l))
    if abs(lz) < value(B)
      return true, lz
    else
      return false, nothing
    end
  else
    return false, nothing
  end
end

#TODO: also allow for
# p-adic, q-adic
# all locals
# finite fields
# rel. ext
# ...
function extension_field(f::ZZPolyRingElem, n::String = "_a"; cached::Bool = true, check::Bool = true)
  return number_field(f, n; cached = cached, check = check)
end
function extension_field(f::QQPolyRingElem, n::String = "_a"; cached::Bool = true, check::Bool = true)
  return number_field(f, n; cached = cached, check = check)
end

function extension_field(f::Generic.Poly{<:Generic.RationalFunctionFieldElem{T}}, n::String = "_a";  cached::Bool = true, check::Bool = true) where {T}
  return function_field(f, n; cached = cached)
end

function extension_field(f::Generic.Poly{AbsSimpleNumFieldElem}, n::String = "_a";  cached::Bool = true, check::Bool = true)
  return number_field(f, n; cached = cached)
end

#Hecke.function_field(f::Generic.Poly{<:Generic.RationalFunctionFieldElem{T}}, n::String = "_a";  cached::Bool = true, check::Bool = true) where {T} = function_field(f, n, cached = cached)


@doc raw"""
Finds a Tschirnhausen transformation, ie a polynomial in `Zx` s.th.
`I` evaluated at the (roots in ) `r` does not have repetitions.

  ``|\{ I^s(t(r_1), ..., t(r_n)) | s in T\}| = |T|``
"""
function find_transformation(r, I::SLPoly, T::AbstractVector{PermGroupElem}; RNG::AbstractRNG=Random.default_rng())
  return find_transformation(r, [I^t for t = T], RNG = RNG)
end

function find_transformation(r, I::Vector{<:SLPoly}; RNG::AbstractRNG = Random.default_rng())
  Zx = Hecke.Globals.Zx
  ts = gen(Zx)
  cnt = 0
  while true
    rt = map(ts, r)
    conj = [evaluate(i, rt) for i = I]
    if length(Set(conj)) == length(I)
      return ts
    end
    while true
      cnt += 1
      cnt > 20 && error("no Tschirni found")
      ts = rand(RNG, Zx, 2:rand(2:max(2, length(r))), -4:4) #TODO: try smaller degrees stronger
      if degree(ts) > 0
        break
      end
    end
  end
end

function relative_invariant(G, U; Chain::Union{Nothing, <:Vector{<:Tuple{PermGroup, <:Any}}} = nothing)
  @assert degree(G) == degree(U)
  S, g = slpoly_ring(ZZ, degree(G), cached = false)
  if index(G, U) == 1 # not type stable
    return one(S), [one(G)]
  end

  if Chain !== nothing
    c = [x[1] for x = Chain]
    I = [x[2] for x = Chain]
    pop!(I)
  else
    c = reverse(maximal_subgroup_chain(G, U))
    I = [invariant(c[i], c[i+1]) for i=1:length(c)-1]
  end

  @vprint :GaloisGroup 2 "using a subgroup chain with orders $(map(order, c))\n"
  #the I[i] is a relative primitive element - it may be absolute or not...
  # right_transversal: U*g, thus G>U>V -> V a b where a runs (U/V), b (G/U):
  # G = cup U b, U = cup V a, so G = V a b
  ts = [gen(Hecke.Globals.Zx) for i = I]
  tv  = [right_transversal(c[i], c[i+1]) for i=1:length(c)-1]
  k = GF(next_prime(10000))
  r = [rand(k) for i=1:degree(G)]
  while length(Set(r)) < length(r)
    r = [rand(k) for i=1:degree(G)]
  end
  mu = ones(Int, length(I))
  #need tschirni per invar
  local conj
  ts[1] = find_transformation(r, I[1], tv[1])
  T = collect(tv[1])
  if ts[1] == gen(Hecke.Globals.Zx)
    a = I[1]
  else
    a = evaluate(I[1], map(ts[1], gens(parent(I[1]))))
  end

  for j=2:length(I)
    local nc, rt
    ts[j] = find_transformation(r, I[j], tv[j])
    lT = length(T)
    for k=2:length(tv[j])
      append!(T, [tv[j][k] * t for t =  T[1:lT]])
    end

    mu[j] = 0

    if ts[j] == gen(Hecke.Globals.Zx)
      b = I[j]
    else
      b = evaluate(I[j], map(ts[j], gens(parent(I[1]))))
    end
    while true
      d = Set([evaluate((mu[j]*a+b)^t, r) for t = T])
      if length(d) == length(T)
        break
      end
      r = [rand(k) for i=1:degree(G)]
      mu[j] += 1
      if mu[j] > characteristic(k)
        error("dnw")
      end
    end
    if iszero(mu[j])
      a = b
    else
      a = mu[j]*a+b
    end
  end
  @vprint :GaloisGroup 2  "have primitive element via $mu \n$a\n"
  return a, T
end

@doc raw"""
    fixed_field(GC::GaloisCtx, U::PermGroup, extra::Int = 5)

Given the `GaloisCtx` as returned by a call to `galois_group` and a subgroup
`U` of the Galois group, compute the field fixed by `U` as a simple
extension.
"""
function fixed_field(GC::GaloisCtx, U::PermGroup, extra::Int = 5)
  G = GC.G
  if index(G, U) == 1
    return Hecke.rationals_as_number_field()[1]
  end
  #XXX: seems to be broken for reducible f, ie. intransitive groups
  a, T = relative_invariant(G, U)
  r = roots(GC, bound_to_precision(GC, GC.B))
  ts = find_transformation(r, a, T, RNG = MersenneTwister(1))

  B = upper_bound(GC, a, ts)
  m = length(T)
  B = m*B^m

  @vtime :GaloisGroup 2 r = roots(GC, bound_to_precision(GC, B, extra))
  if ts != gen(parent(ts))
    r = map(ts, r)
  end
  compile!(a)
  @vtime :GaloisGroup 2 conj = [evaluate(a, t, r) for t = T]

  fl, val = isinteger(GC, B, sum(conj))
  @assert fl
  ps = [val]
  d = copy(conj)
  while length(ps) < m
#    @show length(ps)
    @vtime :GaloisGroup 2 d .*= conj
    fl, val = isinteger(GC, B, sum(d))
    @assert fl
    push!(ps, val)
  end
  ps = map(base_ring(GC.f), ps)

  f = Hecke.power_sums_to_polynomial(ps)
  @assert degree(f) > 0
#  @assert is_irreducible(f)
  k = extension_field(f, check = false, cached = false)[1]
  @assert all(x->isone(denominator(x)), coefficients(k.pol))
  @assert is_monic(k.pol)
  return k
end

"""
    galois_quotient(C::GaloisCtx, Q::PermGroup)

Finds all(?) subfields of the splitting field s.th. the galois
group will be permutation isomorphic to Q.
"""
function galois_quotient(C::GaloisCtx, Q::PermGroup)
  G = C.G
  if order(G) % degree(Q) != 0
    return []
  end
  s = map(representative, subgroup_classes(G, order = divexact(order(G), degree(Q))))
  res = []
  for U = s
    phi = right_coset_action(G, U)
    if is_isomorphic(Q, image(phi)[1])[1]
      push!(res, fixed_field(C, U))
    end
  end
  return res
end

@doc raw"""
    galois_quotient(C::GaloisCtx, d::Int)

Finds all(?) subfields (up to isomorphism) of the splitting field of degree d
with galois group isomorphic to the original one.

# Examples
```jldoctest; filter = r"Galois context\(.*\]\)"
julia> Qx, x = QQ["x"];

julia> G, C = galois_group(x^3-2);

julia> galois_quotient(C, 6)
1-element Vector{Any}:
 Number field of degree 6 over QQ

julia> galois_group(ans[1])
(Permutation group of degree 6 and order 6, Galois context for x^6 + 324*x^4 - 4*x^3 + 34992*x^2 + 1296*x + 1259716 and prime 13)

julia> is_isomorphic(ans[1], G)
true
```
"""
function galois_quotient(C::GaloisCtx, d::Int)
  G = C.G
  if order(G) % d != 0
    return []
  end
  s = map(representative, subgroup_classes(G, order = divexact(order(G), d)))
  res = []
  for U = s
    phi = right_coset_action(G, U)
    if order(image(phi)[1]) == order(G)
      push!(res, fixed_field(C, U))
    end
  end
  return res
end

"""
    galois_quotient(C::GaloisCtx, d::Int, n::Int)

Finds all subfields of the splitting field with galois group the n-th
transitive group in degree d
"""
galois_quotient(C::GaloisCtx, d::Int, n::Int) = 
            galois_quotient(C, transitive_group(d, n))

"""
    galois_quotient(f::PolyRingElem, p::Vector{Int})

Equivalent to

    galois_quotient(galois_group(f)[2], p[1], p[2])

Finds all subfields of the splitting field of f with galois group
the p[2]-th transitive group of degree p[1]
"""
function galois_quotient(f::PolyRingElem, p::Vector{Int})
  G, C = galois_group(f)
  @assert length(p) == 2
  return galois_quotient(C, p[1], p[2])
end

#based on 
#  doi:10.1006/jsco.2000.0376
# Using Galois Ideals for Computing Relative Resolvents
# Aubrey & Valibouze
# JSC (2000) 30, 635-651

function cauchy_ideal(f::ZZPolyRingElem; parent::MPolyRing = polynomial_ring(QQ, degree(f), cached = false)[1])
  return cauchy_ideal(f(gen(Hecke.Globals.Qx)), parent=parent)
end

@doc raw"""
    cauchy_ideal(f::PolyRingElem{<:FieldElem})

The coefficients of `f` are the elementary symmetric functions evaluated
at the roots of `f`. The `cauchy_ideal` is the ideal generated
by the differences between the elementary symmetric functions and the
coefficients.

# Examples
```jldoctest
julia> Qx, x = QQ["x"];

julia> cauchy_ideal(x^4-2)
Ideal generated by
  x4^4 - 2
  x3^3 + x3^2*x4 + x3*x4^2 + x4^3
  x2^2 + x2*x3 + x2*x4 + x3^2 + x3*x4 + x4^2
  x1 + x2 + x3 + x4
```
"""
function cauchy_ideal(f::PolyRingElem{<:FieldElem}; parent::MPolyRing = polynomial_ring(base_ring(f), degree(f), cached = false)[1])
  x = gens(parent)
  n = degree(f)
  f = f(x[n])
  c = [f]
  for i=1:n-1
    f = divexact(evaluate(f, vcat(x[1:n-i], [x[n-i]], x[n-i+2:n])) - f, x[n-i]-x[n-i+1])
    push!(c, f)
  end
  return ideal(c)
end

@doc raw"""
    galois_ideal(C::GaloisCtx, extra::Int = 5)

The so-called Galois ideal is a description of the splitting field of the polynomial
underlying `C`as
a quotient by some maximal ideal. Algebraically, this ideal is an irreducible
component of the Cauchy ideal, the ideal generated by the elementary symmetric
functions and the coefficients of the polynomial.

# Examples
```jldoctest
julia> Qx, x = QQ["x"];

julia> i = galois_ideal(galois_group(x^4-2)[2])
Ideal generated by
  x4^4 - 2
  x3^3 + x3^2*x4 + x3*x4^2 + x4^3
  x2^2 + x2*x3 + x2*x4 + x3^2 + x3*x4 + x4^2
  x1 + x2 + x3 + x4
  x1*x3 + x2*x4
  x1^2*x3^2 + x2^2*x4^2 - 4
  x1^4 - 2
  x2^4 - 2
  x3^4 - 2
  x4^4 - 2

julia> k, _ = number_field(i);

julia> length(roots(k, x^4-2))
4
```
"""
function galois_ideal(C::GaloisCtx, extra::Int = 5)
  #TODO if group is small, then factoring will be better...
  f = C.f
  id = gens(cauchy_ideal(f))
  R = parent(id[1])
  x = gens(R)
  n = length(x)
  #we need to go down (by hand) to the starting group.
  #This is either C.G (if there was no descent in Stauduhar)
  #or the 1st group in the descent chain (C.chn)...
  #TODO: the subfields use, implicitly, special invariants, so
  #      we should be able to avoid the chain
  G = symmetric_group(n)
  if isdefined(C, :start) && C.start[1] == 1 # start with intersection of wreath products
    _, g = slpoly_ring(ZZ, n)
    for bs = C.start[2]
      W = PermGroup(wreath_product(symmetric_group(length(bs[1])), symmetric_group(length(bs))))
      W = W^symmetric_group(n)(vcat(bs...))
      G = intersect(G, W)[1]
      #each subfield causes, possibly, several invariants...
      # (prod(g[b]) for b = bs)
      # are the conjugates of a primitive element, need possibly a shift
      # prod(g[b] .+ i) (according to Klueners)
      # so the elem. symm (or the power sums) of the above need to be in Z
      # furthermore, all roots of this poly are in K, so beta = g(alpha)
      # this gives more invariants... 
      # all of them together basically prove the subfield, so they should
      # be enough for the galois_ideal
      # TODO: if the subfield has a smaller group, then this group
      # could also be incorporated here...(use the galois_ideal of the 
      # subfield and evaluate at the PE below)
      # in general, if G is large, then there are few descents, hence this
      # might work well
      # if G is small, then iterated factoring (possibly using the SolveRadical
      # stuff) is better. See galois_factor there
      r = roots(C, 5)
      k = 0
      while true
        pe = [prod(r[b] .+ k) for b = bs]
        if length(Set(pe)) == length(bs)
          break
        end
        k += 1
      end
      PE = [prod(g[b] .+ k) for b = bs]
      B = upper_bound(C, power_sum, PE, length(PE))
      rt = roots(C, bound_to_precision(C, B))
      pe = [evaluate(I, rt) for I = PE]
      po = pe
      fl, v = isinteger(C, B, sum(pe))
      push!(id, sum(evaluate(I, x) for I = PE) -v)
      @assert fl
      h = [QQ(v)]
      while length(h) < length(PE)
        po = po .* pe
        fl, v = isinteger(C, B, sum(po))
        @assert fl
        push!(h, QQ(v))
        push!(id, sum(evaluate(I^length(h), x) for I = PE) - v)
      end
      h = Hecke.power_sums_to_polynomial(h)
      q = roots(number_field(C.f)[1], h)
      @assert length(q) > 0
      q = parent(defining_polynomial(parent(q[1])))(q[1])
      #TODO: think h(q(y)) is probably boring, while q(y) == pe might be 
      #      not....
      for y = x
        push!(id, h(q(y)))
      end
    end
  end

  if length(C.chn) == 0
    c = maximal_subgroup_chain(G, C.G)
  else
    c = maximal_subgroup_chain(G, C.chn[1][1])
  end

  r = roots(C, bound_to_precision(C, C.B))
  k, mk = residue_field(parent(r[1]))
  r = map(mk, r)

  for i=1:length(c)-1
    I = invariant(c[i+1], c[i]) 
    T = right_transversal(c[i+1], c[i])
    ts = find_transformation(r, I, T)
    B = upper_bound(C, I, ts)
    r = roots(C, bound_to_precision(C, B))
    r = map(ts, r)
    compile!(I)
    for t = T
      e = evaluate(I, t, r)
      fl, v = isinteger(C, B, e)
      if fl
        push!(id, v-evaluate(I, t, map(ts, x)))
        break
      end
    end
  end

  for (_, I, ts, T) = C.chn
    B = upper_bound(C, I, ts)
    r = roots(C, bound_to_precision(C, B, extra))
    r = map(ts, r)
    compile!(I)
    for t = T
      e = evaluate(I, t, r)
      fl, v = isinteger(C, B, e)
      @assert fl
      push!(id, v-evaluate(I, t, x))
    end
  end
  return ideal(id)
end

function Hecke.absolute_primitive_element(K::Oscar.NfNSGen{QQFieldElem, QQMPolyRingElem})
  while true
    a = rand(K, -10:10)
    f = minpoly(a)
    if degree(f) == degree(K)
      return a
    end
  end
end

function Hecke.absolute_minpoly(a::Oscar.NfNSGenElem{QQFieldElem, QQMPolyRingElem})
  return minpoly(a)
end

function blow_up(G::PermGroup, C::GaloisCtx, lf::Vector, con::PermGroupElem=one(G))
  
  if all(x->x[2] == 1, lf)
    return G, C
  end
  
  n = degree(G)
  st = 0
  mp = Dict{Int, Int}(i => i for i=1:n)

  icon = inv(con)

  gs = map(Vector, gens(G))
  for (g, k) = lf
    for j=2:k
      for i=1:degree(g)
        mp[n+i] = (st+i)^con
      end
      for h = gs
        for i=1:degree(g)
          push!(h, h[(st+i)^con]^icon-st+n)
        end
      end
      n += degree(g)
    end
    st += degree(g)
  end

  C.rt_num = mp
  S = symmetric_group(n)
  GG, _ = sub(S, map(S, gs))

  h = hom(G, GG, gens(G), gens(GG))
  @assert is_injective(h) && is_surjective(h)
  return GG, C
end


#TODO: do not move from ZZPolyRingElem to QQPolyRingElem to ZZPolyRingElem...
function galois_group(f::ZZPolyRingElem; pStart::Int = 2*degree(f), prime::Int = 0)
  return galois_group(f(gen(Hecke.Globals.Qx)), pStart = pStart, prime = prime)
end

@doc raw"""
    are_disjoint(G::GaloisCtx, S::GaloisCtx)

Tests if the splitting fields implicitly defined by G and S are
linearly disjoint, or equivalently, if the Galois group for the product
of the polynomials is the direct product of the groups.

If the function returns true, the fields are provable disjoint,
false does not imply no-disjoint.
"""
function are_disjoint(G::GaloisCtx, S::GaloisCtx)
  return false
end
function are_disjoint(G::GaloisCtx{T}, S::GaloisCtx{T}) where T <: Union{Hecke.qAdicRootCtx, ComplexRootCtx}
  #assumes the polynomials (G.f, S.f) are irreducible
  @vprint :GaloisGroup 2 "testing for disjointness...."
  r = isone(gcd(order(G.G), order(S.G)))
  r && @vprint :GaloisGroup 2 "group orders are coprime\n"
  r && return r 
  #if the field discs are coprime, the (splitting) fields are disjoint
  #so if the poly discs are already coprime we are save
  #if the gcd of the poly discs is small, we can try to compute
  #the p-maximal order orders so see if the p-maximal discs are coprime
  #
  
  g =  gcd(ZZ(discriminant(G.f)), ZZ(discriminant(S.f)))
  r = isone(g)
  r && @vprint :GaloisGroup 2 "poly discs are coprime\n"
  if !r #&& (degree(G.f) > 10 || degree(S.f) > 10) && nbits(g) < 20  #try harder
    if degree(G.f) > degree(S.f)
      G, S = S, G
    end 
    @vprint :GaloisGroup 2 "\n poly discs have small gcd ($g)..."
    o1 = any_order(number_field(G.f, cached = false, check = false)[1])
    p = radical(g)
    O1 = pmaximal_overorder(o1, p)
    p = gcd(discriminant(O1), p)
    if isone(p)
      @vprint :GaloisGroup 2 "p-maximal order of 1st is p-free\n"
      return true
    end
    o2 = any_order(number_field(S.f, cached = false, check = false)[1])
    O2 = pmaximal_overorder(o2, p)
    @vprint :GaloisGroup 2 "p-maximal order of 2nd is "
    p = gcd(discriminant(O2), p)
    if isone(p)
      @vprint :GaloisGroup 2 "p-free\n"
      return true
    end
    @vprint :GaloisGroup 2 "not p-free ($p is left), no decision\n"
  end
  return r
end

@doc raw"""
    galois_group(f::PolyRingElem{<:FieldElem})

Computes the automorphism group of a splitting field of `f` as an explicit
group of permutations of the roots. Furthermore, the `GaloisCtx` is
returned allowing algorithmic access to the splitting field.
"""
function galois_group(f::PolyRingElem{<:FieldElem}; prime=0, pStart::Int = 2*degree(f))
  lf = [(k,v) for  (k,v) = factor(f).fac]

  if length(lf) == 1
    @vprint :GaloisGroup 1 "poly has only one factor\n"
    G, C = galois_group(extension_field(lf[1][1], cached = false)[1], pStart = pStart, prime = prime)
    return blow_up(G, C, lf)
  end

  @vprint :GaloisGroup 1 "factoring input...gives\n$lf\n"

  lg = [x[1] for x = lf]
  
  g = prod(lg)
  p, ct = find_prime(g, pStart = pStart, prime = prime)

  C = [galois_group(extension_field(x, cached = false)[1], prime = p)[2] for x = lg]
  #= Strategy:
    sort factors into clusters where the splitting fields are disjoint.
    There the resulting group is the direct product

    Within each cluster the group is a subdirect product
    Maybe use Gap's SubDirectProducts for pairs? It is
    supposed to be better than just doing everything

    Can I can MaximalSubDirectProducts?

    Plan:
     - function to check if 2 GaloisCtx are disjoint. True if disjoiunt,
           false if undecided
     - use this to find the clusters
     - deal with them
     - form the big product
  =#

  g = Graph{Undirected}(length(C))
  for i=1:length(C)
    for j=i+1:length(C)
      if !are_disjoint(C[i], C[j])
        add_edge!(g, i, j)
      end
    end
  end
  cl = connected_components(g)
  @vprint :GaloisGroup 1 "found $(length(cl)) connected components\n"

  res = Vector{Tuple{typeof(C[1]), PermGroupElem}}()
  function setup(C::Vector{<:GaloisCtx})
    G, emb, pro = inner_direct_product([x.G for x = C], morphisms = true)
    g = prod(x.f for x = C)

    CC = GaloisCtx(g, p)
    rr = roots(CC, 5, raw = true) #raw is necessary for non-monic case
                                  #the scaling factor is the lead coeff
                                  #thus not the same for all factors...
    @assert length(Set(rr)) == length(rr)

    d = map(frobenius, rr)
    si = [findfirst(y->y==x, rr) for x = d]

    @vprint :GaloisGroup 1 "found Frobenius element: $si\n"

    k, mk = residue_field(parent(rr[1]))
    rr = map(mk, rr)
    po = Int[]
    for GC = C
      r = roots(GC, 5, raw = true)
      K, mK = residue_field(parent(r[1]))
      r = map(mK, r)
      phi = Hecke.find_morphism(K, k)
      po = vcat(po, [findfirst(x->x == phi(y), rr) for y = r])
    end

    con = (symmetric_group(length(po))(po))
    G = G^con
    CC.G = G #needed as a fallback if no descent is used
    F = GroupFilter()
  #=
    function fi(x)
      @show [pro[y](x^inv(con))[1] == C[y].G for y=1:length(C)]
      return all(y->pro[y](x^inv(con))[1] == C[y].G, 1:length(C))
    end
    push!(F, fi)
  =#  

    push!(F, x->!is_transitive(x), "subdirect case: group is transitive")
    push!(F, x->all(y->pro[y](x^inv(con))[1] == C[y].G, 1:length(C)), "subdirect case: wrong projections")
    return CC, G, F, G(si), con
  end
  for X = cl
    @vprint :GaloisGroup 1 "dealing with factors $X ...\n"
    if length(X) == 1
      push!(res, (C[X][1], one(C[X][1].G)))
    else
      function red(A::GaloisCtx, B::GaloisCtx)
        CC, G, F, fr, con = setup([A, B])
        return descent(CC, G, F, fr, grp_id = x->(:-, :-))[2], con
      end
      local con
      #=
      G = C2 x C2 x C2 x C2 e.g. sqrt(2) sqrt(3) sqrt(5) sqrt(30)
           (a,a)     (a,a)     with a: group combinations that have to be
          a         a                  testet
          a              a     with b: groups that are not sub-direct enough
              a     a
              a          a
          a   a     a
          a   a          a
          a   a     a    a
          a         a    a
              a     a    a

          b
             b
                    b
                         b
          b  b
                    b    b
      the idea is tp combine pairs in the hope that some descents can be 
      found on smaller groups (smaller support/ degree)
      We'll have to test all possible subdirect subgroups, but this tries 
      smaller cases first
      =#
      @vprint :GaloisGroup 2 "entering combination tree..."
      while length(X) > 1
        local p
        d = map(x->degree(x.f), C[X])
        p = sortperm(d)
        @vprint :GaloisGroup 2 "combining $(X[p[1:2]]) of degree $(degree(C[X[p[1]]].f)) and $(degree(C[X[p[2]]].f))\n"
        Hecke.pushindent()
        D, con = red(C[X[p[1]]], C[X[p[2]]])
        Hecke.popindent()
        #delete the larger index as to not upset the smaller one
        #then replace the Ctx at the smaller index
        if p[1] < p[2]
          deleteat!(X, p[2])
          C[X[p[1]]] = D
        else
          deleteat!(X, p[1])
          C[X[p[2]]] = D
        end
      end
      push!(res, (C[X[1]], con))
    end
  end

  @vprint :GaloisGroup 1 "combining clusters...\n"
  if length(res) == 1
    CC = res[1][1]
    G = CC.G
    con = res[1][2]
  else
    CC, G, F, fr, con = setup([x[1] for x = res])
  end
  @vprint :GaloisGroup 1 "adding multiplicity...\n"

  return blow_up(G, CC, lf, con)
end

function Nemo.cyclotomic(n::Int, x::QQPolyRingElem)
  return Nemo.cyclotomic(n, gen(Hecke.Globals.Zx))(x)
end


################################################################################
#
#  Promote rules
#
################################################################################

AbstractAlgebra.promote_rule(::Type{BoundRingElem}, ::Type{ZZRingElem}) = BoundRingElem

AbstractAlgebra.promote_rule(::Type{BoundRingElem}, ::Type{T}) where {T <: Integer} = BoundRingElem

include("Group.jl")
include("POSet.jl")
include("SeriesEval.jl")
include("Qt.jl")
include("RelGalois.jl")

end

using .GaloisGrp

#=
       M12: 2-transitive, hence msum is a waste
x^12 - 4*x^11 + 4*x^10 + 12*x^9 - 72*x^8 + 168*x^7 - 132*x^6 - 324*x^5 + 1197*x^4 - 1752*x^3 + 1500*x^2 - 672*x + 207 

several index 2 groups, uses generic invar
x^16 + 209*x^14 + 102*x^13 + 18138*x^12 + 13232*x^11 + 814855*x^10 + 673869*x^9 + 19699456*x^8 + 17373605*x^7 + 258261711*x^6 + 283233089*x^5 + 2368948579*x^4 + 2075376015*x^3 + 13009666150*x^2 + 14279830429*x + 37222748001


primitive, 2-transitive
x^10 - 2*x^9 - 263*x^8 + 1136*x^7 + 18636*x^6 - 120264*x^5 - 81916*x^4 + 1314656*x^3 - 290197*x^2 - 3135542*x + 2052019

needs the group theory test for isinvar in line 805 and related
x^16 + 4*x^15 + 20*x^14 + 44*x^13 + 106*x^12 + 120*x^11 + 180*x^10 + 44*x^9 + 134*x^8 - 120*x^7 + 172*x^6 - 32*x^5 + 598*x^4 + 312*x^3 + 616*x^2 + 147
=#