oscar-system / Oscar.jl
 1 #module MPolyModule  2 3 import AbstractAlgebra: PolyRing, PolynomialRing, total_degree, degree, Ideal,  4  MPolyElem, Generic.MPolyExponentVectors, Generic.MPolyCoeffs,  5  Generic.MPolyBuildCtx, Generic.push_term!, Generic.finish, MPolyRing,  6  base_ring, ngens, gens, dim, ordering, SetMap, Map  7 import Nemo  8 import Nemo: fmpz, fmpq  9 10 import Singular  11 12 import Hecke  13 import Hecke: MapHeader, math_html  14 15 16 17 export PolynomialRing, total_degree, degree, MPolyIdeal, MPolyElem, ordering, ideal, coordinates,  18  jacobi_matrix, jacobi_ideal, normalize, divrem, isprimary, isprime  19 20 ##############################################################################  21 #  22 # could/ should be in AbstractAlgebra  23 #  24 # some sugar to make creation of strange rings easier  25 # possibly lacks the passing of the ordering...  26 ##############################################################################  27 28 #TODO: reduce = divrem in Nemo. Should be faster - if we have the correct basis  29 30 #allows  31 # PolynomialRing(QQ, :a=>1:3, "b"=>1:3, "c=>1:5:10)  32 # -> QQx, [a1, a2, a3], [b1 ,b2, b3], ....  33 10 function PolynomialRing(R::AbstractAlgebra.Ring, v::Pair{<:Union{String, Symbol}, <:Union{StepRange{Int, Int}, UnitRange{Int}}}...; cached::Bool = false, ordering::Symbol = :lex)  34 12  s = String[]  35 12  g = []  36 10  j = 1  37 12  for (a, b) = v  38 12  h = []  39 12  for i = b  40 12  if occursin('#', "$a")  41 0  aa = replace("$a", '#' => "$i")  42  else  43 12  if Hecke.inNotebook()  44 0  aa = "$(a)_{$i}"  45  else  46 12  aa = "$a$i"  47  end  48  end  49 12  push!(s, aa)  50 12  push!(h, j)  51 12  j += 1  52  end  53 12  push!(g, h)  54  end  55 12  Rx, c = PolynomialRing(R, s, cached = cached, ordering = ordering)  56 12  return Rx, [c[x] for x = g]...  57 end  58 59 10 function Base.getindex(R::MPolyRing, i::Int)  60 12  i == 0 && return zero(R)  61 12  return gen(R, i)  62 end  63 64 ######################################################################  65 # pretty printing for iJulia notebooks..  66 #  67 68 0 function Base.show(io::IO, mime::IJuliaMime, R::MPolyRing)  69 0  io = IOContext(io, :compact => true)  70 0  print(io, "\$")  71 0  math_html(io, R)  72 0  print(io, "\$")  73 end  74 75 0 function math_html(io::IO, R::MPolyRing)  76 0  print(io, "\\text{Multivariate Polynomial Ring in$(nvars(R)) variables:} ")  77 0  math_html(io, gens(R))  78 0  print(io, "\\text{ over }")  79 0  math_html(io, base_ring(R))  80 end  81 82 0 function math_html(io::IO, R::MPolyElem)  83 0  f = "$R"  84 0  f = replace(f, r"_\$([0-9]*)" => s"t_{\1}")  85 0  f = replace(f, "*" => "")  86 0  f = replace(f, r"\^([0-9]*)" => s"^{\1}")  87 0  print(io, f)  88 end  89 90 0 function Base.show(io::IO, ::IJuliaMime, R::MPolyElem)  91 0  print(io, "\$")  92 0  math_html(io, R)  93 0  print(io, "\$")  94 end  95 96 ###################################################  97 98 module Orderings  99 100 using Oscar, Markdown  101 import Oscar: Ring, MPolyRing, MPolyElem, weights  102 export anti_diagonal, lex, degrevlex, deglex, weights, MonomialOrdering, singular  103 104 abstract type AbsOrdering end  105 """  106 Ring-free monomial ordering: just the indices of the variables are given.  107 T can be a UnitRange to make Singular happy or any Array if the  108  variables are not consequtive  109 """  110 mutable struct GenOrdering{T} <: AbsOrdering  111  vars::T  112  ord::Symbol  113  wgt::fmpz_mat  114 0  function GenOrdering(u::T, s::Symbol) where {T <: AbstractArray{Int, 1}}  115 0  r = new{typeof(u)}()  116 0  r.vars = u  117 0  r.ord = s  118 0  return r  119  end  120 0  function GenOrdering(u::T, m::fmpz_mat; ord::Symbol = :weight) where {T <: AbstractArray{Int, 1}}  121 0  r = new{typeof(u)}()  122 0  @assert ncols(m) == length(u)  123 0  r.vars = u  124 0  r.ord = ord  125 0  r.wgt = m  126 0  return r  127  end  128 end  129 130 """  131 The product of a and b (vcat of the the matrices)  132 """  133 mutable struct ProdOrdering <: AbsOrdering  134  a::AbsOrdering  135  b::AbsOrdering  136 end  137 138 0 Base.:*(a::AbsOrdering, b::AbsOrdering) = ProdOrdering(a, b)  139 140 #not really user facing  141 0 function ordering(a::AbstractArray{Int, 1}, s::Union{Symbol, fmpz_mat})  142 0  i = minimum(a)  143 0  I = maximum(a)  144 0  if I-i+1 == length(a) #testif variables are consecutive or not.  145 0  return GenOrdering(i:I, s)  146  end  147 0  return GenOrdering(collect(a), s)  148 end  149 150 #not really user facing  151 0 function ordering(a::AbstractArray{Int, 1}, s::Symbol, w::fmpz_mat)  152 0  i = minimum(a)  153 0  I = maximum(a)  154 0  if I-i+1 == length(a)  155 0  return GenOrdering(i:I, w, ord = s)  156  end  157 0  return GenOrdering(collect(a), w, ord = s)  158 end  159 160 161 #not really user facing, flattens a product of product orderings into an array  162 0 function flat(a::GenOrdering)  163 0  return [a]  164 end  165 0 function flat(a::ProdOrdering)  166 0  return vcat(flat(a.a), flat(a.b))  167 end  168 169 @doc Markdown.doc"""  170  anti_diagonal(R::Ring, n::Int)  171 172 A square matrix with 1 on the anti-diagonal.  173 """  174 0 function anti_diagonal(R::Ring, n::Int)  175 0  a = zero_matrix(R, n, n)  176 0  for i=1:n  177 0  a[i, n-i+1] = one(R)  178  end  179 0  return a  180 end  181 182 #not user facing  183 0 function weights(a::GenOrdering)  184 0  if a.ord == :lex || a.ord == Symbol("Singular(lp)")  185 0  return identity_matrix(ZZ, length(a.vars))  186  end  187 0  if a.ord == :deglex  188 0  return [matrix(ZZ, 1, length(a.vars), ones(fmpz, length(a.vars)));  189  identity_matrix(ZZ, length(a.vars)-1) zero_matrix(ZZ, length(a.vars)-1, 1)]  190  end  191 0  if a.ord == :degrevlex || a.ord == Symbol("Singular(dp)")  192 0  return [matrix(ZZ, 1, length(a.vars), ones(fmpz, length(a.vars))) ;  193  zero_matrix(ZZ, length(a.vars)-1, 1) anti_diagonal(ZZ, length(a.vars)-1)]  194  end  195 0  if a.ord == Symbol("Singular(ls)")  196 0  return -identity_matrix(ZZ, length(a.vars))  197  end  198 0  if a.ord == Symbol("Singular(ds)")  199 0  return [-matrix(ZZ, 1, length(a.vars), ones(fmpz, length(a.vars))) ;  200  zero_matrix(ZZ, length(a.vars)-1, 1) anti_diagonal(ZZ, length(a.vars)-1)]  201  end  202 0  if a.ord == Symbol("Singular(a)") || a.ord == Symbol("Singular(M)")  203 0  return a.wgt  204  end  205 end  206 207 #not user facing  208 0 function weights(a::AbsOrdering)  209 0  aa = flat(a)  210 0  m = matrix(ZZ, 0, 0, [])  211 0  for o = aa  212 0  w = weights(o)  213 0  if maximum(o.vars) > ncols(m)  214 0  m = hcat(m, zero_matrix(ZZ, nrows(m), maximum(o.vars) - ncols(m)))  215  end  216 0  mm = zero_matrix(ZZ, nrows(w), ncols(m))  217 0  for r = 1:nrows(w)  218 0  for c = 1:length(o.vars)  219 0  mm[r, o.vars[c]] = w[r, c]  220  end  221  end  222 0  m = vcat(m, mm)  223  end  224 0  return m  225 end  226 227 """  228 Orderings actually applied to polynomial rings (as opposed to variable indices)  229 """  230 mutable struct MonomialOrdering{S}  231  R::S  232  o::AbsOrdering  233 end  234 235 #not really user facing, not exported  236 @doc Markdown.doc"""  237  ordering(a::Vector{MPolyElem}, s::Symbol)  238  ordering(a::Vector{MPolyElem}, m::fmpz_mat)  239  ordering(a::Vector{MPolyElem}, s::Symbol, m::fmpz_mat)  240 241 Defines an ordering to be applied to the variables in a.  242 In the first form the symbol s has to be one of :lex, :deglex or :degrevlex.  243 In the second form, a weight ordering using the given matrix is used.  244 In the last version, the symbol if of the form Singular(..).  245 """  246 0 function ordering(a::AbstractArray{<:MPolyElem, 1}, s...)  247 0  R = parent(first(a))  248 0  g = gens(R)  249 0  aa = [findfirst(x -> x == y, g) for y = a]  250 0  if nothing in aa  251 0  error("only variables allowed")  252  end  253 0  return ordering(aa, s...)  254 end  255 256 @doc Markdown.doc"""  257  :*(M::MonomialOrdering, N::MonomialOrdering)  258 259 For orderings on the same ring, the product ordering obained by concatenation  260 of the weight matrics.  261 """  262 0 function Base.:*(M::MonomialOrdering, N::MonomialOrdering)  263 0  M.R == N.R || error("wrong rings")  264 0  return MonomialOrdering(M.R, M.o*N.o)  265 end  266 267 0 function Base.show(io::IO, M::MonomialOrdering)  268 0  a = flat(M.o)  269 0  if length(a) > 1  270 0  print(io, "Product ordering: ")  271 0  for i=1:length(a)-1  272 0  show(io, M.R, a[i])  273 0  print(io, " \\times ")  274  end  275  end  276 0  show(io, M.R, a[end])  277 end  278 279 0 function Base.show(io::IO, R::MPolyRing, o::GenOrdering)  280 0  if o.ord == :weight  281 0  print(io, "weight($(gens(R)[o.vars]) via$(o.wgt))")  282  else  283 0  print(io, "$(String(o.ord))($(gens(R)[o.vars]))")  284  end  285 end  286 287 @doc Markdown.doc"""  288  lex(v::AbstractArray{<:MPolyElem, 1}) -> MonomialOrdering  289 290 Defines the lex (lexicographic) ordering on the variables given.  291 """  292 0 function lex(v::AbstractArray{<:MPolyElem, 1})  293 0  return MonomialOrdering(parent(first(v)), ordering(v, :lex))  294 end  295 @doc Markdown.doc"""  296  deglex(v::AbstractArray{<:MPolyElem, 1}) -> MonomialOrdering  297 298 Defines the deglex ordering on the variables given.  299 """  300 0 function deglex(v::AbstractArray{<:MPolyElem, 1})  301 0  return MonomialOrdering(parent(first(v)), ordering(v, :deglex))  302 end  303 @doc Markdown.doc"""  304  degrevlex(v::AbstractArray{<:MPolyElem, 1}) -> MonomialOrdering  305 306 Defines the degreveex ordering on the variables given.  307 """  308 0 function degrevlex(v::AbstractArray{<:MPolyElem, 1})  309 0  return MonomialOrdering(parent(first(v)), ordering(v, :degrevlex))  310 end  311 312 @doc Markdown.doc"""  313  singular(ord::Symbol, v::AbstractArray{<:MPolyElem, 1}) -> MonomialOrdering  314 315 Defines an ordering given in terms of Singular primitives on the variables given.  316 ord can be one of :lp, :ls, :dp, :ds.  317 """  318 0 function singular(ord::Symbol, v::AbstractArray{<:MPolyElem, 1})  319 0  return MonomialOrdering(parent(first(v)), ordering(v, Symbol("Singular($(string(ord)))")))  320 end  321 322 @doc Markdown.doc"""  323  singular(ord::Symbol, v::AbstractArray{<:MPolyElem, 1}, w::AbstractArray{Int, 2}) -> MonomialOrdering  324 325 Defines an ordering given in terms of Singular weight ordering (M) with the  326 matrix given. ord has to be :M here.  327 """  328 0 function singular(ord::Symbol, v::AbstractArray{<:MPolyElem, 1}, w::Array{<:Union{Integer, fmpz}, 2})  329 0  @assert ord == :M  330 0  W = matrix(ZZ, size(w, 1), size(w, 2), w)  331 0  return MonomialOrdering(parent(first(v)), ordering(v, Symbol("Singular($(string(ord)))"), W))  332 end  333 334 @doc Markdown.doc"""  335  singular(ord::Symbol, v::AbstractArray{<:MPolyElem, 1}, w::AbstractArray{Int, 2}) -> MonomialOrdering  336 337 Defines an ordering given in terms of Singular weight ordering (a) with the  338 weights given. ord has to be :a here. The weights will be supplemented by  339 0.  340 """  341 0 function singular(ord::Symbol, v::AbstractArray{<:MPolyElem, 1}, w::Array{<:Union{Integer, fmpz}, 1})  342 0  @assert ord == :a  343 0  W = map(fmpz, w)  344 0  while length(v) > length(W)  345 0  push!(W, 0)  346  end  347 348 0  return MonomialOrdering(parent(first(v)), ordering(v, Symbol("Singular($(string(ord)))"), matrix(ZZ, 1, length(W), W)))  349 350 end  351 352 @doc Markdown.doc"""  353  weights(M::MonomialOrdering)  354   355 Compute a corresponding weight matrix for the given ordering.  356 """  357 0 function weights(M::MonomialOrdering)  358 0  return weights(M.o)  359 end  360 361 @doc Markdown.doc"""  362  simplify(M::MonomialOrdering) -> MonomialOrdering  363 364 Compute a weight ordering with a unique weight matrix.  365 """  366 0 function simplify(M::MonomialOrdering)  367 0  w = weights(M)  368 0  ww = matrix(ZZ, 0, ncols(w), [])  369 0  for i=1:nrows(w)  370 0  if iszero_row(w, i)  371 0  continue  372  end  373 0  nw = w[i, :]  374 0  c = content(nw)  375 0  if c != 1  376 0  nw = divexact(nw, c)  377  end  378 0  for j=1:nrows(ww)  379 0  h = findfirst(x->ww[j, x] != 0, 1:ncols(w))  380 0  if nw[1, h] != 0  381 0  nw = abs(ww[j, h])*nw - sign(ww[j, h])*nw[1, h]*ww[j, :]  382  end  383  end  384 0  if !iszero(nw)  385 0  c = content(nw)  386 0  if !isone(c)  387 0  nw = divexact(nw, c)  388  end  389 0  ww = vcat(ww, nw)  390  end  391  end  392 0  return MonomialOrdering(M.R, ordering(1:ncols(ww), ww))  393 end  394 395 import Base.==  396 0 function ==(M::MonomialOrdering, N::MonomialOrdering)  397 0  return simplify(M).o.wgt == simplify(N).o.wgt  398 end  399 400 0 function Base.hash(M::MonomialOrdering, u::UInt)  401 0  return hash(simplify(M).o.wgt, u)  402 end  403 404 end # module Orderings  405 406 using .Orderings  407 export lex, deglex, degrevlex, weights, MonomialOrdering, singular  408 409 410 411 ##############################################################################  412 #  413 # workhorse: BiPolyArray  414 # ideals are (mostly) generated on the Nemo side, but structural computations  415 # are in Singular. To avoid permanent conversion, the list of generators = sideal  416 # is captured in BiPolyArray: for Ocsar this is Array{MPoly, 1}  417 # Singular sideal  418 #  419 #TODO/ to think  420 # default in Nemo is :lex  421 # Singular is :degrevlex -> better for std  422 #by default, use different orders???  423 #make BiPolyArray use different orders for both? make the type depend on it?  424 #  425 #for std: abstraction to allow Christian to be used  426 #  427 #type for orderings, use this...  428 #in general: all algos here needs revision: do they benefit from gb or not?  429 430 mutable struct BiPolyArray{S}  431  O::Array{S, 1}  432  S::Singular.sideal  433  Ox #Oscar Poly Ring  434  Sx # Singular Poly Ring, poss. with different ordering  435  isGB::Bool #if the Singular side (the sideal) will be a GB  436  ord :: Orderings.AbsOrdering #for this ordering  437 10  function BiPolyArray(a::Array{T, 1}; keep_ordering::Bool = true, isGB::Bool = false) where {T <: MPolyElem}  438 12  r = new{T}()  439 12  r.O = a  440 12  r.Ox = parent(a[1])  441 12  r.isGB = isGB  442 12  r.Sx = singular_ring(r.Ox, keep_ordering = keep_ordering)  443 12  return r  444  end  445 10  function BiPolyArray(Ox::T, b::Singular.sideal) where {T <: MPolyRing}  446 12  r = new{elem_type(T)}()  447 12  r.S = b  448 12  r.O = Array{elem_type(T)}(undef, Singular.ngens(b))  449 12  r.Ox = Ox  450 12  r.isGB = b.isGB  451 12  r.Sx = base_ring(b)  452 12  return r  453  end  454 end  455 456 0 function Base.getindex(A::BiPolyArray, ::Val{:S}, i::Int)  457 0  if !isdefined(A, :S)  458 0  A.S = Singular.Ideal(A.Sx, [A.Sx(x) for x = A.O])  459  end  460 0  return A.S[i]  461 end  462 463 10 function Base.getindex(A::BiPolyArray, ::Val{:O}, i::Int)  464 12  if !isassigned(A.O, i)  465 12  A.O[i] = A.Ox(A.S[i])  466  end  467 12  return A.O[i]  468 end  469 470 function Base.length(A::BiPolyArray)  471 12  if isdefined(A, :S)  472 12  return Singular.ngens(A.S)  473  else  474 12  return length(A.O)  475  end  476 end  477 478 function Base.iterate(A::BiPolyArray, s::Int = 1)  479 12  if s > length(A)  480 12  return nothing  481  end  482 12  return A[Val(:O), s], s+1  483 end  484 485 10 Base.eltype(::BiPolyArray{S}) where S = S  486 487 ##############################################################################  488 #  489 # Conversion to and from Singular: in particular, some Rings are  490 # special as they exist natively in Singular and thus should be used  491 #  492 ##############################################################################  493 #  494 # Needs convert(Target(Ring), elem)  495 # Ring(s.th.)  496 #  497 # singular_ring(Nemo-Ring) tries to create the appropriate Ring  498 #  499 500 10 function (Ox::MPolyRing)(f::Singular.spoly)  501 12  O = base_ring(Ox)  502 12  Sx = parent(f)  503 12  @assert ngens(Sx) == ngens(Ox)  504 12  g = MPolyBuildCtx(Ox)  505 12  for (c, e) = Base.Iterators.zip(Singular.coefficients(f), Singular.exponent_vectors(f))  506 12  push_term!(g, O(c), e)  507  end  508 12  return finish(g)  509 end  510 511 512 10 function (S::Singular.Rationals)(a::fmpq)  513 12  b = Base.Rational{BigInt}(a)  514 12  return S(b)  515 end  516 517 12 (F::Singular.N_ZpField)(a::Nemo.gfp_elem) = F(lift(a))  518 12 (F::Singular.N_ZpField)(a::Nemo.nmod) = F(lift(a))  519 12 (F::Nemo.GaloisField)(a::Singular.n_Zp) = F(Int(a))  520 12 (F::Nemo.NmodRing)(a::Singular.n_Zp) = F(Int(a))  521 522 #Note: Singular crashes if it gets Nemo.ZZ instead of Singular.ZZ ((Coeffs(17)) instead of (ZZ))  523 12 singular_ring(::Nemo.FlintIntegerRing) = Singular.Integers()  524 12 singular_ring(::Nemo.FlintRationalField) = Singular.Rationals()  525 12 singular_ring(F::Nemo.GaloisField) = Singular.Fp(Int(characteristic(F)))  526 12 singular_ring(F::Nemo.NmodRing) = Singular.Fp(Int(characteristic(F)))  527 0 singular_ring(R::Singular.PolyRing; keep_ordering::Bool = true) = R  528 529 10 function singular_ring(Rx::MPolyRing{T}; keep_ordering::Bool = true) where {T <: RingElem}  530 12  if keep_ordering  531 12  return Singular.PolynomialRing(singular_ring(base_ring(Rx)),  532  [string(x) for x = Nemo.symbols(Rx)],  533  ordering = ordering(Rx),  534  cached = false)[1]  535  else  536 12  return Singular.PolynomialRing(singular_ring(base_ring(Rx)),  537  [string(x) for x = Nemo.symbols(Rx)],  538  cached = false)[1]  539  end  540 end  541 542 10 function singular_ring(Rx::MPolyRing{T}, ord::Symbol) where {T <: RingElem}  543 12  return Singular.PolynomialRing(singular_ring(base_ring(Rx)),  544  [string(x) for x = Nemo.symbols(Rx)],  545  ordering = ord,  546  cached = false)[1]  547 end  548 549 #catch all for generic nemo rings  550 0 function Oscar.singular_ring(F::AbstractAlgebra.Ring)  551 0  return Singular.CoefficientRing(F)  552 end  553 554 0 function (b::AbstractAlgebra.Ring)(a::Singular.n_unknown)  555 0  Singular.libSingular.julia(Singular.libSingular.cast_number_to_void(a.ptr))  556 end  557 558 ##############################################################################  559 #  560 # Multivariate ideals - also used for the decorated stuff  561 #  562 ##############################################################################  563 564 mutable struct MPolyIdeal{S} <: Ideal{S}  565  gens::BiPolyArray{S}  566  gb::BiPolyArray{S}  567  dim::Int  568 569  function MPolyIdeal(g::Array{T, 1}) where {T <: MPolyElem}  570 12  r = new{T}()  571 12  r.dim = -1 #not known  572 12  r.gens = BiPolyArray(g, keep_ordering = false)  573 12  return r  574  end  575 10  function MPolyIdeal(Ox::T, s::Singular.sideal) where {T <: MPolyRing}  576 12  r = new{elem_type(T)}()  577 12  r.dim = -1 #not known  578 12  r.gens = BiPolyArray(Ox, s)  579 12  if s.isGB  580 0  r.gb = r.gens  581  end  582 12  return r  583  end  584  function MPolyIdeal(B::BiPolyArray{T}) where T  585 0  r = new{T}()  586 0  r.dim = -1  587 0  r.gens = B  588 0  return r  589  end  590 end  591 592 4 function Base.show(io::IO, I::MPolyIdeal)  593 4  print(io, "ideal generated by: ")  594 4  g = collect(I.gens)  595 4  first = true  596 4  for i = g  597 4  if first  598 4  print(io, i)  599 4  first = false  600  else  601 4  print(io, ", ", i)  602  end  603  end  604 4  print(io, "")  605 end  606 607 0 function Base.show(io::IO, ::IJuliaMime, I::MPolyIdeal)  608 0  print(io, "\$")  609 0  math_html(io, I)  610 0  print(io, "\$")  611 end  612 613 0 function math_html(io::IO, I::MPolyIdeal)  614 0  print(io, "\\text{ideal generated by: }")  615 0  g = collect(I.gens)  616 0  first = true  617 0  for i = g  618 0  if first  619 0  math_html(io, i)  620 0  first = false  621  else  622 0  print(io, ", ")  623 0  math_html(io, i)  624  end  625  end  626 0  print(io, "")  627 end  628 629 630 631 632 633 10 function ideal(g::Array{Any, 1})  634 12  return ideal(typeof(g[1])[x for x = g])  635 end  636 637 638 639 10 function ideal(Rx::MPolyRing, s::Singular.sideal)  640 12  return MPolyIdeal(Rx, s)  641 end  642 643 function singular_assure(I::MPolyIdeal)  644 12  singular_assure(I.gens)  645 end  646 647 10 function singular_assure(I::BiPolyArray)  648 12  if !isdefined(I, :S)  649 12  I.S = Singular.Ideal(I.Sx, [I.Sx(x) for x = I.O])  650 12  if I.isGB  651 12  I.S.isGB = true  652  end  653  end  654 end  655 656 657 0 function oscar_assure(I::MPolyIdeal)  658 0  if !isdefined(I.gens, :O)  659 0  I.gens.O = [I.gens.Ox(x) for x = gens(I.gens.S)]  660  end  661 end  662 663 0 function Base.copy(f::MPolyElem)  664 0  Ox = parent(f)  665 0  g = MPolyBuildCtx(Ox)  666 0  for (c,e) = Base.Iterators.zip(MPolyCoeffs(f), MPolyExponentVectors(f))  667 0  push_term!(g, c, e)  668  end  669 0  return finish(g)  670 end  671 672 10 function map_entries(R, M::Singular.smatrix)  673 12  s = nrows(M), ncols(M)  674 12  S = parent(R(zero(base_ring(M))))  675 12  return matrix(S, s[1], s[2], elem_type(S)[R(M[i,j]) for i=1:s[1] for j=1:s[2]])  676 end  677 678 0 function syzygy_module(a::Array{MPolyElem, 1})  679  #only graded modules exist  680 0  error("not implemented yet")  681 end  682 683 4 function (F::Generic.FreeModule)(s::Singular.svector)  684 4  pv = Tuple{Int, elem_type(base_ring(F))}[]  685 4  pos = Int[]  686 4  values = []  687 4  Rx = base_ring(F)  688 4  R = base_ring(Rx)  689 4  for (i, e, c) = s  690 4  f = Base.findfirst(x->x==i, pos)  691 4  if f === nothing  692 4  push!(values, MPolyBuildCtx(Rx))  693 4  f = length(values)  694 4  push!(pos, i)  695  end  696 4  push_term!(values[f], R(c), e)  697  end  698 4  pv = Tuple{Int, elem_type(Rx)}[(pos[i], finish(values[i])) for i=1:length(pos)]  699 4  e = zero(F)  700 4  for (k,v) = pv  701 4  e += v*gen(F, k)  702  end  703 4  return e  704 end  705 706 @doc Markdown.doc"""  707  jacobi_matrix(f::MPolyElem)  708 709 Given a polynomial$f$this function returns the Jacobian matrix J_f=(\partial_{x_1}f,...,\partial_{x_n}f)^T of$f$.  710 """  711 10 function jacobi_matrix(f::MPolyElem)  712 12  R = parent(f)  713 12  n = nvars(R)  714 12  return matrix(R, n, 1, [derivative(f, i) for i=1:n])  715 end  716 717 @doc Markdown.doc"""  718  jacobi_ideal(f::MPolyElem)  719 720 Given a polynomial$f$this function returns the Jacobian ideal of$f$.  721 """  722 10 function jacobi_ideal(f::MPolyElem)  723 12  R = parent(f)  724 12  n = nvars(R)  725 12  return ideal(R, [derivative(f, i) for i=1:n])  726 end  727 728 @doc Markdown.doc"""  729  jacobi_matrix(g::Array{<:MPolyElem, 1})  730 731 Given an array g=[f_1,...,f_m] of polynomials over the same base ring,  732 this function returns the Jacobian matrix J=(\partial_{x_i}f_j)_{i,j} of g.  733 """  734 10 function jacobi_matrix(g::Array{<:MPolyElem, 1})  735 12  R = parent(g[1])  736 12  n = nvars(R)  737 12  @assert all(x->parent(x) == R, g)  738 12  return matrix(R, n, length(g), [derivative(x, i) for i=1:n for x = g])  739 end  740 741 ##########################################  742 #  743 # Singular library related functions  744 #  745 ##########################################  746 747 ##########################  748 #  749 # basic maps  750 #  751 ##########################  752 0 function im_func(f::MPolyElem, S::MPolyRing, i::Array{Int, 1})  753 0  O = base_ring(S)  754 0  g = MPolyBuildCtx(S)  755 0  for (c, e) = Base.Iterators.zip(MPolyCoeffs(f), MPolyExponentVectors(f))  756 0  f = zeros(Int, nvars(S))  757 0  for j=1:length(e)  758 0  if i[j] == 0  759 0  e[j] != 0 && error("illegal map: var$(j) is used")  760  else  761 0  f[i[j]] = e[j]  762  end  763  end  764 0  push_term!(g, O(c), f)  765  end  766 0  return finish(g)  767 end  768 769 770 abstract type OscarMap <: SetMap end  771 772 mutable struct MPolyHom_vars{T1, T2} <: Map{T1, T2, Hecke.HeckeMap, MPolyHom_vars}  773  header::Hecke.MapHeader  774  Hecke.@declare_other  775  i::Array{Int, 1}  776 777 0  function MPolyHom_vars{T1, T2}(R::T1, S::T2, i::Array{Int, 1}) where {T1 <: MPolyRing, T2 <: MPolyRing}  778 0  r = new()  779 0  p = sortperm(i)  780 0  j = Int[]  781 0  for h = 1:length(p)  782 0  if i[p[h]] != 0  783 0  j = p[h:length(p)]  784 0  break  785  end  786  end  787 0  r.header = MapHeader{T1, T2}(R, S, x -> im_func(x, S, i), y-> im_func(y, R, j))  788 0  r.i = i  789 0  return r  790  end  791 792 0  function MPolyHom_vars{T1, T2}(R::T1, S::T2; type::Symbol = :none) where {T1 <: MPolyRing, T2 <: MPolyRing}  793 794 0  if type == :names  795 0  i = Int[]  796 0  for h = symbols(R)  797 0  push!(i, findfirst(x -> x == h, symbols(S)))  798  end  799 0  return MPolyHom_vars{T1, T2}(R, S, i)  800  end  801 0  error("type not supported")  802  end  803 end  804 805 0 (f::MPolyHom_vars)(g::MPolyElem) = image(f, g)  806 807 function Hecke.hom(R::MPolyRing, S::MPolyRing, i::Array{Int, 1})  808 0  return MPolyHom_vars{typeof(R), typeof(S)}(R, S, i)  809 end  810 811 0 function _lift(S::Singular.sideal, T::Singular.sideal)  812 0  R = base_ring(S)  813 0  @assert base_ring(T) == R  814 0  c, r = Singular.libSingular.id_Lift(S.ptr, T.ptr, R.ptr)  815 0  M = Singular.Module(R, c)  816 817 0  if Singular.ngens(M) == 0 || iszero(M[1])  818 0  error("elem not in module")  819  end  820 0  return M  821 end  822 823 #TODO: return a matrix??  824 @doc Markdown.doc"""  825  coordinates(a::Array{<:MPolyElem, 1}, b::Array{<:MPolyElem, 1})  826 827 Tries to write the entries of b as linear combinations of a.  828 """  829 0 function coordinates(a::Array{<:MPolyElem, 1}, b::Array{<:MPolyElem, 1})  830 0  ia = ideal(a)  831 0  ib = ideal(b)  832 0  singular_assure(ia)  833 0  singular_assure(ib)  834 0  c = _lift(ia.gens.S, ib.gens.S)  835 0  F = free_module(parent(a[1]), length(a))  836 0  return [F(c[x]) for x = 1:Singular.ngens(c)]  837 end  838 839 0 function coordinates(a::Array{<:MPolyElem, 1}, b::MPolyElem)  840 0  return coordinates(a, [b])[1]  841 end  842 843 ###################################################  844 845 # Some isless functions for orderings:  846 # _isless_:ord(f, k, l) returns true if the k-th term is lower than the l-th  847 # term of f in the ordering :ord.  848 849 10 function _isless_lex(f::MPolyElem, k::Int, l::Int)  850 12  n = nvars(parent(f))  851 12  for i = 1:n  852 12  ek = exponent(f, k, i)  853 12  el = exponent(f, l, i)  854 12  if ek == el  855 0  continue  856 12  elseif ek > el  857 12  return false  858  else  859 12  return true  860  end  861  end  862 0  return false  863 end  864 865 10 function _isless_neglex(f::MPolyElem, k::Int, l::Int)  866 12  n = nvars(parent(f))  867 12  for i = 1:n  868 12  ek = exponent(f, k, i)  869 12  el = exponent(f, l, i)  870 12  if ek == el  871 0  continue  872 12  elseif ek < el  873 0  return false  874  else  875 12  return true  876  end  877  end  878 0  return false  879 end  880 881 10 function _isless_revlex(f::MPolyElem, k::Int, l::Int)  882 12  n = nvars(parent(f))  883 12  for i = n:-1:1  884 12  ek = exponent(f, k, i)  885 12  el = exponent(f, l, i)  886 12  if ek == el  887 0  continue  888 12  elseif ek > el  889 0  return false  890  else  891 12  return true  892  end  893  end  894 0  return false  895 end  896 897 10 function _isless_negrevlex(f::MPolyElem, k::Int, l::Int)  898 12  n = nvars(parent(f))  899 12  for i = n:-1:1  900 12  ek = exponent(f, k, i)  901 12  el = exponent(f, l, i)  902 12  if ek == el  903 0  continue  904 12  elseif ek < el  905 12  return false  906  else  907 0  return true  908  end  909  end  910 0  return false  911 end  912 913 10 function _isless_deglex(f::MPolyElem, k::Int, l::Int)  914 12  tdk = total_degree(term(f, k))  915 12  tdl = total_degree(term(f, l))  916 12  if tdk < tdl  917 12  return true  918 0  elseif tdk > tdl  919 0  return false  920  end  921 0  return _isless_lex(f, k, l)  922 end  923 924 10 function _isless_degrevlex(f::MPolyElem, k::Int, l::Int)  925 12  tdk = total_degree(term(f, k))  926 12  tdl = total_degree(term(f, l))  927 12  if tdk < tdl  928 12  return true  929 0  elseif tdk > tdl  930 0  return false  931  end  932 0  return _isless_negrevlex(f, k, l)  933 end  934 935 10 function _isless_negdeglex(f::MPolyElem, k::Int, l::Int)  936 12  tdk = total_degree(term(f, k))  937 12  tdl = total_degree(term(f, l))  938 12  if tdk > tdl  939 0  return true  940 12  elseif tdk < tdl  941 12  return false  942  end  943 0  return _isless_lex(f, k, l)  944 end  945 946 10 function _isless_negdegrevlex(f::MPolyElem, k::Int, l::Int)  947 12  tdk = total_degree(term(f, k))  948 12  tdl = total_degree(term(f, l))  949 12  if tdk > tdl  950 0  return true  951 12  elseif tdk < tdl  952 12  return false  953  end  954 0  return _isless_negrevlex(f, k, l)  955 end  956 957 # Returns the degree of the k-th term of f weighted by w,  958 # that is deg(x^a) = w_1a_1 + ... + w_na_n.  959 # No sanity checks are performed!  960 function weighted_degree(f::MPolyElem, k::Int, w::Vector{Int})  961 12  ek = exponent_vector(f, k)  962 12  return dot(ek, w)  963 end  964 965 10 function _isless_weightlex(f::MPolyElem, k::Int, l::Int, w::Vector{Int})  966 12  dk = weighted_degree(f, k, w)  967 12  dl = weighted_degree(f, l, w)  968 12  if dk < dl  969 0  return true  970 12  elseif dk > dl  971 0  return false  972  end  973 12  return _isless_lex(f, k, l)  974 end  975 976 10 function _isless_weightrevlex(f::MPolyElem, k::Int, l::Int, w::Vector{Int})  977 12  dk = weighted_degree(f, k, w)  978 12  dl = weighted_degree(f, l, w)  979 12  if dk < dl  980 0  return true  981 12  elseif dk > dl  982 0  return false  983  end  984 12  return _isless_negrevlex(f, k, l)  985 end  986 987 0 function _isless_weightneglex(f::MPolyElem, k::Int, l::Int, w::Vector{Int})  988 0  dk = weighted_degree(f, k, w)  989 0  dl = weighted_degree(f, l, w)  990 0  if dk < dl  991 0  return true  992 0  elseif dk > dl  993 0  return false  994  end  995 0  return _isless_lex(f, k, l)  996 end  997 998 0 function _isless_weightnegrevlex(f::MPolyElem, k::Int, l::Int, w::Vector{Int})  999 0  dk = weighted_degree(f, k, w)  1000 0  dl = weighted_degree(f, l, w)  1001 0  if dk > dl  1002 0  return true  1003 0  elseif dk < dl  1004 0  return false  1005  end  1006 0  return _isless_negrevlex(f, k, l)  1007 end  1008 1009 10 function _isless_matrix(f::MPolyElem, k::Int, l::Int, M::Union{ Array{T, 2}, MatElem{T} }) where T  1010 12  ek = exponent_vector(f, k)  1011 12  el = exponent_vector(f, l)  1012 12  n = nvars(parent(f))  1013 12  for i = 1:size(M, 1)  1014 12  eki = sum( M[i, j]*ek[j] for j = 1:n )  1015 12  eli = sum( M[i, j]*el[j] for j = 1:n )  1016 12  if eki == eli  1017 0  continue  1018 12  elseif eki > eli  1019 0  return false  1020  else  1021 12  return true  1022  end  1023  end  1024 0  return false  1025 end  1026 1027 10 function _perm_of_terms(f::MPolyElem, ord_lt::Function)  1028 12  p = collect(1:length(f))  1029 12  sort!(p, lt = (k, l) -> ord_lt(f, k, l), rev = true)  1030 12  return p  1031 end  1032 1033 # Requiring R for consistence with the other lt_from_ordering functions  1034 function lt_from_ordering(::MPolyRing, ord::Symbol)  1035 12  if ord == :lex || ord == :lp  1036 12  return _isless_lex  1037 12  elseif ord == :revlex || ord == :rp  1038 12  return _isless_revlex  1039 12  elseif ord == :deglex || ord == :Dp  1040 12  return _isless_deglex  1041 12  elseif ord == :degrevlex || ord == :dp  1042 12  return _isless_degrevlex  1043 12  elseif ord == :neglex || ord == :ls  1044 12  return _isless_neglex  1045 12  elseif ord == :negrevlex || ord == :rs  1046 12  return _isless_negrevlex  1047 12  elseif ord == :negdeglex || ord == :Ds  1048 12  return _isless_negdeglex  1049 12  elseif ord == :negdegrevlex || ord == :ds  1050 12  return _isless_negdegrevlex  1051  else  1052 0  error("Ordering $ord not available")  1053  end  1054 end  1055 1056 10 function lt_from_ordering(R::MPolyRing, ord::Symbol, w::Vector{Int})  1057 12  @assert length(w) == nvars(R) "Number of weights has to match number of variables"  1058 1059 12  if ord == :weightlex || ord == :Wp  1060 12  @assert all(x -> x > 0, w) "Weights have to be positive"  1061 12  return (f, k, l) -> _isless_weightlex(f, k, l, w)  1062 12  elseif ord == :weightrevlex || ord == :wp  1063 12  @assert all(x -> x > 0, w) "Weights have to be positive"  1064 12  return (f, k, l) -> _isless_weightrevlex(f, k, l, w)  1065 0  elseif ord == :weightneglex || ord == :Ws  1066 0  @assert !iszero(w[1]) "First weight must not be 0"  1067 0  return (f, k, l) -> _isless_weightneglex(f, k, l, w)  1068 0  elseif ord == :weightnegrevlex || ord == :ws  1069 0  @assert !iszero(w[1]) "First weight must not be 0"  1070 0  return (f, k, l) -> _isless_weightnegrevlex(f, k, l, w)  1071  else  1072 0  error("Ordering$ord not available")  1073  end  1074 end  1075 1076 function lt_from_ordering(R::MPolyRing, M::Union{ Array{T, 2}, MatElem{T} }) where T  1077 12  @assert size(M, 2) == nvars(R) "Matrix dimensions have to match number of variables"  1078 1079 12  return (f, k, l) -> _isless_matrix(f, k, l, M)  1080 end  1081 1082 10 function terms(f::MPolyElem, ord::Function)  1083 12  perm = _perm_of_terms(f, ord)  1084 12  return ( term(f, perm[i]) for i = 1:length(f) )  1085 end  1086 1087 10 function coefficients(f::MPolyElem, ord::Function)  1088 12  perm = _perm_of_terms(f, ord)  1089 12  return ( coeff(f, perm[i]) for i = 1:length(f) )  1090 end  1091 1092 10 function exponent_vectors(f::MPolyElem, ord::Function)  1093 12  perm = _perm_of_terms(f, ord)  1094 12  return ( exponent_vector(f, perm[i]) for i = 1:length(f) )  1095 end  1096 1097 10 function monomials(f::MPolyElem, ord::Function)  1098 12  perm = _perm_of_terms(f, ord)  1099 12  return ( monomial(f, perm[i]) for i = 1:length(f) )  1100 end  1101 1102 for s in (:terms, :coefficients, :exponent_vectors, :monomials)  1103  @eval begin  1104 10  function ($s)(f::MPolyElem, ord::Symbol)  1105 12  R = parent(f)  1106 12  if ord == ordering(R)  1107 12  return ($s)(f)  1108  end  1109 1110 12  lt = lt_from_ordering(R, ord)  1111 12  return ($s)(f, lt)  1112  end  1113 1114 10  function ($s)(f::MPolyElem, M::Union{ Array{T, 2}, MatElem{T} }) where T  1115 12  R = parent(f)  1116 12  lt = lt_from_ordering(R, M)  1117 12  return ($s)(f, lt)  1118  end  1119 1120 10  function ($s)(f::MPolyElem, ord::Symbol, weights::Vector{Int})  1121 12  R = parent(f)  1122 12  lt = lt_from_ordering(R, ord, weights)  1123 12  return ($s)(f, lt)  1124  end  1125 1126 0  function ($s)(f::MPolyElem, ord::MonomialOrdering)  1127 0  R = parent(f)  1128 0  lt = lt_from_ordering(R, weights(ord))  1129 0  return ($s)(f, lt)  1130  end  1131  end  1132 end  1133 1134 for s in ("term", "coefficient", "monomial")  1135  @eval begin  1136 10  function ($(Symbol("leading_$s")))(args...)  1137 12  return first($(Symbol("\$(s)s"))(args...))  1138  end  1139  end  1140 end  1141 1142 0 function leading_term(f::MPolyElem)  1143 0  return leading_term(f, ordering(parent(f)))  1144 end  1145 1146 0 function leading_coefficient(f::MPolyElem)  1147 0  return leading_coefficient(f, ordering(parent(f)))  1148 end  1149 1150 0 function leading_monomial(f::MPolyElem)  1151 0  return leading_monomial(f, ordering(parent(f)))  1152 end  1153 1154 ##############################################################################  1155 #  1156 ##############################################################################  1157 1158 #=  1159 function factor(f::MPolyElem)  1160  I = ideal(parent(f), [f])  1161  fS = Singular.factor(I.gens[Val(:S), 1])  1162  R = parent(f)  1163  return Nemo.Fac(R(fS.unit), Dict(R(k) =>v for (k,v) = fS.fac))  1164 end  1165 =#  1166 1167 ################################################################################  1168 1169 @doc Markdown.doc"""  1170  divrem(a::Vector{T}, b::Vector{T}) where T <: MPolyElem{S} where S <: RingElem  1171 Return an array of tuples (qi, ri) consisting of an array of polynomials qi, one  1172 for each polynomial in b, and a polynomial ri such that  1173 a[i] = sum_i b[i]*qi + ri.  1174 """  1175 10 function divrem(a::Vector{T}, b::Vector{T}) where T <: MPolyElem{S} where S <: RingElem  1176 12  return [divrem(x, b) for x in a]  1177 end  1178 1179 ################################################################################  1180 1181 10 function Base.:*(f::MPolyElem, I::MPolyIdeal)  1182 12  R = base_ring(I)  1183 12  R == parent(f) || error("base rings do not match")  1184 12  return ideal(R, f.*gens(I))  1185 end  1186 1187 ################################################################################  1188

Read our documentation on viewing source code .