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 .

Loading