Allow field coercion and pre-set coefficient fields for polyhedra

Project.toml

@@ -36,7 +36,7 @@ AlgebraicSolving = "0.10.1"
 Compat = "4.13.0"
 Distributed = "1.6"
 GAP = "0.16.2"
-Hecke = "0.39.4"
+Hecke = "0.39.7"
 JSON = "1.0.1"
 JSON3 = "1.13.2"
 LazyArtifacts = "1.6"

docs/src/PolyhedralGeometry/intro.md

@@ -86,6 +86,14 @@ coefficient_field(x::PolyhedralObject)
 embedded_number_field
 ```
 
+If there is a corresponding conversion for the field elements, then polyhedra can
+be converted to a different field, e.g., a polytope over a number field can be
+converted to a polytope over the algebraic closure of the rationals.
+
+```@docs
+polyhedron(f::scalar_type_or_field, p::Polyhedron)
+```
+
 ## Type compatibility
 
 When working in polyhedral geometry it can prove advantageous to have various

src/PolyhedralGeometry/PolyhedralGeometry.jl

@@ -1,4 +1,4 @@
-const AnyVecOrMat = Union{MatElem,AbstractVecOrMat}
+const AnyVecOrMat = Union{MatElem,SMat,AbstractVecOrMat}
 
 include("helpers.jl")
 include("iterators.jl")

src/PolyhedralGeometry/Polyhedron/constructors.jl

@@ -31,24 +31,53 @@ polyhedron(A; non_redundant::Bool=false, is_bounded::Union{Bool,Nothing}=nothing
   )
 
 @doc raw"""
-    polyhedron(P::Polymake.BigObject)
+    polyhedron([f::scalar_type_or_field,] P::Polymake.BigObject)
 
 Construct a `Polyhedron` corresponding to a `Polymake.BigObject` of type
-`Polytope`. Scalar type and parent field will be detected automatically. To
+`Polytope`.
+By default, scalar type and parent field will be detected automatically. To
 improve type stability and performance, please use
 [`Polyhedron{T}(p::Polymake.BigObject, f::Field) where T<:scalar_types`](@ref)
 instead, where possible.
+Passing an explicit field will force a conversion (or copy) of all elements to
+that field.
 """
 function polyhedron(p::Polymake.BigObject)
   T, f = _detect_scalar_and_field(Polyhedron, p)
   if T == EmbeddedNumFieldElem{AbsSimpleNumFieldElem} &&
     Hecke.is_quadratic_type(number_field(f))[1] &&
     Polymake.bigobject_eltype(p) == "QuadraticExtension"
-    p = _polyhedron_qe_to_on(p, f)
+    p = _polyhedron_coerce_field(p, f)
   end
   return Polyhedron{T}(p, f)
 end
 
+function polyhedron(f::scalar_type_or_field, p::Polymake.BigObject)
+  T = f isa Field ? elem_type(f) : f
+  p = _polyhedron_coerce_field(p, f)
+  return Polyhedron{T}(p, f)
+end
+
+@doc raw"""
+    polyhedron(f::scalar_type_or_field, P::Polyhedron)
+
+Try to coerce a given polyhedron `P` to a different field `f`.
+
+# Examples
+```jldoctest
+julia> d = dodecahedron()
+Polytope in ambient dimension 3 with EmbeddedAbsSimpleNumFieldElem type coefficients
+
+julia> dqb = polyhedron(algebraic_closure(QQ), d)
+Polytope in ambient dimension 3 with QQBarFieldElem type coefficients
+```
+"""
+function polyhedron(f::scalar_type_or_field, p::Polyhedron)
+  f, T = f isa Field ? (f, elem_type(f)) : _determine_parent_and_scalar(f)
+  p = _polyhedron_coerce_field(pm_object(p), f)
+  return Polyhedron{T}(p, f)
+end
+
 @doc raw"""
     polyhedron([::Union{Type{T}, Field},] A::AnyVecOrMat, b) where T<:scalar_types
 

src/PolyhedralGeometry/Polyhedron/standard_constructions.jl

@@ -253,25 +253,28 @@ end
 cube(d::Int, l, u) = cube(QQFieldElem, d, l, u)
 
 @doc raw"""
-    tetrahedron()
+    tetrahedron([F::scalar_type_or_field])
 
 Construct the regular tetrahedron, one of the Platonic solids.
 """
 tetrahedron() = polyhedron(Polymake.polytope.tetrahedron());
+tetrahedron(f::scalar_type_or_field) = polyhedron(F, Polymake.polytope.tetrahedron());
 
 @doc raw"""
-    dodecahedron()
+    dodecahedron([F::scalar_type_or_field])
 
 Construct the regular dodecahedron, one out of two Platonic solids.
 """
 dodecahedron() = polyhedron(Polymake.polytope.dodecahedron());
+dodecahedron(F::scalar_type_or_field) = polyhedron(F, Polymake.polytope.dodecahedron());
 
 @doc raw"""
-    icosahedron()
+    icosahedron([F::scalar_type_or_field])
 
 Construct the regular icosahedron, one out of two exceptional Platonic solids.
 """
 icosahedron() = polyhedron(Polymake.polytope.icosahedron());
+icosahedron(F::scalar_type_or_field) = polyhedron(F, Polymake.polytope.icosahedron());
 
 const _johnson_indexes_from_oscar = Set{Int}([9, 10, 13, 16, 17, 18, 20, 21, 22, 23, 24,
   25, 30, 32, 33, 34, 35, 36, 38, 39, 40,
@@ -282,26 +285,34 @@ const _johnson_indexes_from_oscar = Set{Int}([9, 10, 13, 16, 17, 18, 20, 21, 22,
   90, 92])
 
 @doc raw"""
-    johnson_solid(i::Int)
+    johnson_solid([F::scalar_type_or_field,] i::Int)
 
-Construct the `i`-th proper Johnson solid.
+Construct the `i`-th proper Johnson solid. Optionally embedded in a given field, if possible.
 
 A Johnson solid is a 3-polytope whose facets are regular polygons, of various gonalities.
 It is proper if it is not an Archimedean solid.  Up to scaling there are exactly 92 proper Johnson solids.
   See the [Polytope Wiki](https://polytope.miraheze.org/wiki/List_of_Johnson_solids)
 
 See also [`is_johnson_solid`](@ref is_johnson_solid(P::Polyhedron)).
 """
-function johnson_solid(index::Int)
+johnson_solid(index::Int) = _johnson_solid(index)
+johnson_solid(F::scalar_type_or_field, index::Int) = _johnson_solid(index, F)
+
+function _johnson_solid(index::Int, F::Union{scalar_type_or_field,Nothing}=nothing)
   if index in _johnson_indexes_from_oscar
     # code used for generation of loaded files can be found at:
     # https://github.com/dmg-lab/RegularSolidsSrc
     str_index = lpad(index, 2, '0')
     filename = "j$str_index" * ".mrdi"
-    return load(joinpath(oscardir, "data", "JohnsonSolids", filename))
+    pmp = load(joinpath(oscardir, "data", "JohnsonSolids", filename))
+  else
+    pmp = Polymake.polytope.johnson_solid(index)
+    if F === nothing
+      # autodetect field type
+      pmp = polyhedron(pmp)
+    end
   end
-  pmp = Polymake.polytope.johnson_solid(index)
-  return polyhedron(pmp)
+  return F !== nothing ? polyhedron(F, pmp) : pmp
 end
 
 @doc raw"""
@@ -746,7 +757,7 @@ end
 cross_polytope(d::Int64) = cross_polytope(QQFieldElem, d)
 
 @doc raw"""
-    platonic_solid(s)
+    platonic_solid([F::scalar_type_or_field,] s)
 
 Construct a Platonic solid with the name given by String `s` from the list
 below.
@@ -777,14 +788,16 @@ julia> n_facets(T)
 ```
 """
 platonic_solid(s::String) = polyhedron(Polymake.polytope.platonic_solid(s))
+platonic_solid(F::scalar_type_or_field, s::String) =
+  polyhedron(F, Polymake.polytope.platonic_solid(s))
 
 const _archimedean_strings_from_oscar = Set{String}(["snub_cube", "snub_dodecahedron"])
 
 @doc raw"""
-    archimedean_solid(s)
+    archimedean_solid([F::scalar_type_or_field,] s)
 
 Construct an Archimedean solid with the name given by String `s` from the list
-below.
+below. Optionally embedded in a given field, if possible.
 
 See also [`is_archimedean_solid`](@ref is_archimedean_solid(P::Polyhedron)).
 
@@ -838,15 +851,23 @@ julia> n_facets(T)
 14
 ```
 """
-function archimedean_solid(s::String)
+archimedean_solid(s::String) = _archimedean_solid(s)
+archimedean_solid(F::scalar_type_or_field, s::String) = _archimedean_solid(s, F)
+
+function _archimedean_solid(s::String, F::Union{scalar_type_or_field,Nothing}=nothing)
   if s in _archimedean_strings_from_oscar
     filename = s * ".mrdi"
     # code used for generation of loaded files can be found at:
     # https://github.com/dmg-lab/RegularSolidsSrc
-    return load(joinpath(oscardir, "data", "ArchimedeanSolids", filename))
+    pmp = load(joinpath(oscardir, "data", "ArchimedeanSolids", filename))
+  else
+    pmp = Polymake.polytope.archimedean_solid(s)
+    if F === nothing
+      # autodetect field type
+      pmp = polyhedron(pmp)
+    end
   end
-  pmp = Polymake.polytope.archimedean_solid(s)
-  return polyhedron(pmp)
+  return F !== nothing ? polyhedron(F, pmp) : pmp
 end
 
 @doc raw"""
@@ -878,9 +899,10 @@ const _catalan_strings_from_oscar = Set{String}([
 ])
 
 @doc raw"""
-    catalan_solid(s::String)
+    catalan_solid([F::scalar_type_or_field,] s::String)
 
 Construct a Catalan solid with the name `s` from the list below.
+Optionally embedded in a given field, if possible.
 
 # Arguments
 - `s::String`: The name of the desired Archimedean solid.
@@ -931,15 +953,23 @@ julia> n_facets(T)
 12
 ```
 """
-function catalan_solid(s::String)
+catalan_solid(s::String) = _catalan_solid(s)
+catalan_solid(F::scalar_type_or_field, s::String) = _catalan_solid(s, F)
+
+function _catalan_solid(s::String, F::Union{scalar_type_or_field,Nothing}=nothing)
   if s in _catalan_strings_from_oscar
     filename = s * ".mrdi"
     # code used for generation of loaded files can be found at:
     # https://github.com/dmg-lab/RegularSolidsSrc
-    return load(joinpath(oscardir, "data", "CatalanSolids", filename))
+    pmp = load(joinpath(oscardir, "data", "CatalanSolids", filename))
+  else
+    pmp = Polymake.polytope.catalan_solid(s)
+    if F === nothing
+      # autodetect field type
+      pmp = polyhedron(pmp)
+    end
   end
-  pmp = Polymake.polytope.catalan_solid(s)
-  return polyhedron(pmp)
+  return F !== nothing ? polyhedron(F, pmp) : pmp
 end
 
 @doc raw"""

src/PolyhedralGeometry/helpers.jl

@@ -341,6 +341,15 @@ function (NF::Hecke.EmbeddedNumField)(x::Polymake.QuadraticExtension{Polymake.Ra
   return convert(QQFieldElem, g.a) + convert(QQFieldElem, g.b) * a
 end
 
+function (qqb::QQBarField)(x::Polymake.QuadraticExtension{Polymake.Rational})
+  g = Polymake.generating_field_elements(x)
+  if g.r == 0 || g.b == 0
+    return qqb(convert(QQFieldElem, g.a))
+  end
+  r = qqb(convert(QQFieldElem, g.r))
+  return convert(QQFieldElem, g.a) + convert(QQFieldElem, g.b) * sqrt(r)
+end
+
 (F::Field)(x::Polymake.Rational) = F(QQ(x))
 (F::Field)(x::Polymake.OscarNumber) = F(Polymake.unwrap(x))
 
@@ -427,10 +436,12 @@ homogenize(field, mat::AbstractMatrix, val::Union{Number,scalar_types_extended}=
   augment(field, mat, fill(val, size(mat, 1)))
 homogenize(field, mat::MatElem, val::Union{Number,scalar_types_extended}=1) =
   homogenize(field, Matrix(mat), val)
+homogenize(field, mat::SMat, val::Union{Number,scalar_types_extended}=1) =
+  homogenize(field, Matrix(mat), val)
 homogenize(field, nothing, val::Union{Number,scalar_types_extended}) = nothing
 homogenized_matrix(
   field,
-  x::Union{AbstractVecOrMat,MatElem,Nothing},
+  x::Union{AbstractVecOrMat,MatElem,SMat,Nothing},
   val::Union{Number,scalar_types_extended},
 ) =
   homogenize(field, x, val)
@@ -589,11 +600,12 @@ _find_elem_type(x::Any) = typeof(x)
 _find_elem_type(x::Type) = x
 _find_elem_type(x::Polymake.Rational) = QQFieldElem
 _find_elem_type(x::Polymake.Integer) = ZZRingElem
-_find_elem_type(x::AbstractArray) = reshape(_find_elem_type.(x), :)
+_find_elem_type(x::AbstractArray) = collect(reshape(_find_elem_type.(x), :))
 _find_elem_type(x::Tuple) = reduce(vcat, _find_elem_type.(x))
 _find_elem_type(x::AbstractArray{<:AbstractArray}) =
   reduce(vcat, _find_elem_type.(x); init=[])
 _find_elem_type(x::MatElem) = [elem_type(base_ring(x))]
+_find_elem_type(x::SMat) = [elem_type(base_ring(x))]
 
 function _guess_fieldelem_type(x...)
   types = filter(!=(Any), _find_elem_type(x))
@@ -651,6 +663,9 @@ end
 function _determine_parent_and_scalar(::Type{QQFieldElem}, x...)
   return (QQ, QQFieldElem)
 end
+function _determine_parent_and_scalar(::Type{QQBarFieldElem}, x...)
+  return (algebraic_closure(QQ), QQBarFieldElem)
+end
 
 function _determine_parent_and_scalar(::Type{T}, x...) where {T<:scalar_types}
   p = _parent_or_coefficient_field(T, x...)
@@ -825,7 +840,11 @@ function Polymake._fieldelem_to_rational(e::EmbeddedNumFieldElem)
   return Rational{BigInt}(QQ(e))
 end
 
-function Polymake._fieldelem_is_rational(e::EmbeddedNumFieldElem)
+function Polymake._fieldelem_to_rational(e::QQBarFieldElem)
+  return Rational{BigInt}(e)
+end
+
+function Polymake._fieldelem_is_rational(e::Union{EmbeddedNumFieldElem,QQBarFieldElem})
   return is_rational(e)
 end
 
@@ -843,29 +862,44 @@ end
 
 # convert a Polymake.BigObject's scalar from QuadraticExtension to OscarNumber (Polytope only)
 
-function _polyhedron_qe_to_on(x::Polymake.BigObject, f::Field)
-  res = Polymake.polytope.Polytope{Polymake.OscarNumber}()
+function _polyhedron_coerce_field(x::Polymake.BigObject, f::scalar_type_or_field)
+  f = f isa AbstractAlgebra.Floats{Float64} ? Float64 : f
+  T = f isa Field ? elem_type(f) : f
+  res = Polymake.polytope.Polytope{_scalar_type_to_polymake(T)}()
   for pn in Polymake.list_properties(x)
     prop = Polymake.give(x, pn)
-    Polymake.take(res, string(pn), _property_qe_to_on(prop, f))
+    Polymake.take(res, string(pn), _property_coerce_field(prop, f))
   end
   return res
 end
 
-_property_qe_to_on(x::Polymake.BigObject, f::Field) =
+_property_coerce_field(x::Polymake.BigObject, f::scalar_type_or_field) =
   Polymake.BigObject(Polymake.bigobject_type(x), x)
 
-_property_qe_to_on(x::Polymake.PropertyValue, f::Field) = x
+_property_coerce_field(x::Polymake.PropertyValue, f::scalar_type_or_field) = x
 
-_property_qe_to_on(x::Polymake.QuadraticExtension{Polymake.Rational}, f::Field) = f(x)
+_property_coerce_field(
+  x::Polymake.QuadraticExtension{Polymake.Rational}, f::scalar_type_or_field
+) = f(x)
 
-function _property_qe_to_on(x, f::Field)
-  if hasmethod(length, (typeof(x),)) &&
-    eltype(x) <: Polymake.QuadraticExtension{Polymake.Rational}
-    return f.(x)
-  else
-    return x
+_property_coerce_field(x::Polymake.Rational, f::scalar_type_or_field) = f(x)
+
+_property_coerce_field(x::Polymake.OscarNumber, f::scalar_type_or_field) =
+  f(Polymake.unwrap(x))
+
+function _property_coerce_field(x, f::scalar_type_or_field)
+  if hasmethod(length, (typeof(x),))
+    # we need to make sure the output is a proper polymake matrix
+    if eltype(x) <: Polymake.OscarNumber
+      res = similar(x, _scalar_type_to_polymake(f isa Field ? elem_type(f) : f))
+      res .= f.([Polymake.unwrap(e) for e in x])
+      return res
+    elseif eltype(x) <:
+      Union{Polymake.QuadraticExtension{Polymake.Rational},Polymake.Rational}
+      return f.(x)
+    end
   end
+  return x
 end
 
 # Helper function for conversion

test/PolyhedralGeometry/polyhedral_complex.jl

@@ -76,7 +76,7 @@
     @test lineality_space(PCL) isa SubObjectIterator{RayVector{T}}
     @test length(lineality_space(PCL)) == 1
     @test lineality_space(PCL) == [L[:]]
-    @test generator_matrix(lineality_space(PCL)) == matrix(QQ, L)
+    @test generator_matrix(lineality_space(PCL)) == matrix(f, L)
 
     @test lineality_dim(PCFL) == 1
     @test f_vector(PCL) == [0, 4, 4, 1]