Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Better pbw ideals arith #1591

Merged
merged 3 commits into from
Oct 5, 2022
Merged

Better pbw ideals arith #1591

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
d69fa47
pbw_algebra: correct ideal * and ^
tthsqe12 Oct 4, 2022
393408d
check the conditions by default
tthsqe12 Oct 4, 2022
660febb
change error message
tthsqe12 Oct 5, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
131 changes: 70 additions & 61 deletions src/Rings/PBWAlgebra.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@ end

mutable struct PBWAlgIdeal{D, T, S}
basering::PBWAlgRing{T, S}
sdata::Singular.sideal{Singular.spluralg{S}} # the gens of this ideal
sdata::Singular.sideal{Singular.spluralg{S}} # the gens of this ideal, always defined
sopdata::Singular.sideal{Singular.spluralg{S}} # the gens mapped to the opposite
gb::Singular.sideal{Singular.spluralg{S}}
opgb::Singular.sideal{Singular.spluralg{S}}
Expand All @@ -35,8 +35,12 @@ mutable struct PBWAlgIdeal{D, T, S}
d.isTwoSided = (D == 0)
return new{D, T, S}(p, d)
end
function PBWAlgIdeal{D, T, S}(p::PBWAlgRing{T, S}) where {D, T, S}
return new{D, T, S}(p)
function PBWAlgIdeal{D, T, S}(p::PBWAlgRing{T, S},
d::Singular.sideal{Singular.spluralg{S}},
opd::Singular.sideal{Singular.spluralg{S}}) where {D, T, S}
d.isTwoSided = (D == 0)
opd.isTwoSided = (D == 0)
return new{D, T, S}(p, d, opd)
end
end
# the meaning of the direction parameter D
Expand Down Expand Up @@ -292,7 +296,7 @@ end

####
@doc Markdown.doc"""
pbw_algebra(R::MPolyRing{T}, rel, ord::MonomialOrdering) where T
pbw_algebra(R::MPolyRing{T}, rel, ord::MonomialOrdering; check = true) where T

Comment on lines 298 to 300
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is this now checked by default? Isn't this horribly expensive?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Perhaps we should discuss this in general and set a standard: speed over correctness or correctness over speed? Both have clear pros and cons...

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Personally I prefer correct, slow results over wrong, fast results... but of course we all prefer correct & fast results... The question is whether to expect people to explicitly say "trust me, I am an expert" (so it's their own fault for running against a wall; but it our fault if things are "slow" if one does not know about check=false). On the other hand, will people who are more cautious expect our functions to not validate inputs by default? Hmmm...

So yeah, my preference actually is to have check=true the default everywhere and then expect everyone to explicitly say check=false if they are sure about it. But I also realize I may be a minority here :-)

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think that if the default is check = false, it makes the check argument useless. Then the only reason one would want to call a function with check = true is if one expects that the input is not valid. But then one should call the function which validates the input and not wait for an error in the constructor.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am not suggesting to start with check = false as default. Here is another example from the section on modules:

Warning

The functions do not check whether the resulting homomorphism is well-defined, that is, whether it sends the relations of M into the relations of N.

If you are uncertain with regard to well-definedness, use the function below. Note, however, that the check performed by the function requires a Gröbner basis computation. This may take some time.

is_welldefined(a::ModuleFPHom)

Return true if a is well-defined, and false otherwise.

I agree that we should agree on discussing this in general.

In the current example, we could do some time checks first.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Isn't this horribly expensive?

No, it is just O(nvars^3) multiplications and additions in the algebra.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

O.k., thanks!

Given a multivariate polynomial ring `R` over a field, say ``R=K[x_1, \dots, x_n]``, given
a strictly upper triangular matrix `rel` with entries in `R` of type ``c_{ij} \cdot x_ix_j+d_{ij}``,
Expand All @@ -310,7 +314,8 @@ A = K\langle x_1, \dots , x_n \mid x_jx_i = c_{ij} \cdot x_ix_j+d_{ij}, \ 1\leq
See the definition of PBW-algebras in the OSCAR documentation for details.

!!! warning
The `K`-basis condition above is not checked by the function.
The `K`-basis condition above is checked by default. This check may be
skipped by passing `check = false`.

# Examples
```jldoctest
Expand All @@ -324,7 +329,7 @@ julia> A, (x,y,z) = pbw_algebra(R, REL, deglex(gens(R)))
(PBW-algebra over Rational Field in x, y, z with relations y*x = x*y, z*x = x*z, z*y = y*z + 1, PBWAlgElem{fmpq, Singular.n_Q}[x, y, z])
```
"""
function pbw_algebra(r::MPolyRing{T}, rel, ord::MonomialOrdering) where T
function pbw_algebra(r::MPolyRing{T}, rel, ord::MonomialOrdering; check = true) where T
n = nvars(r)
nrows(rel) == n && ncols(rel) == n || error("oops")
scr = singular_coeff_ring(coefficient_ring(r))
Expand All @@ -342,6 +347,9 @@ function pbw_algebra(r::MPolyRing{T}, rel, ord::MonomialOrdering) where T
srel[i,j] = t
end
s, gs = Singular.GAlgebra(sr, C, D)
if check && !iszero(Singular.LibNctools.ndcond(s))
error("PBW-basis condition not satisfied")
end
R = PBWAlgRing{T, S}(s, srel, coefficient_ring(r))
return R, [PBWAlgElem(R, x) for x in gs]
end
Expand All @@ -352,7 +360,7 @@ function weyl_algebra(K::Ring, xs::Vector{Symbol}, dxs::Vector{Symbol})
n == length(dxs) || error("number of differentials should match number of variables")
r, v = PolynomialRing(K, vcat(xs, dxs))
rel = elem_type(r)[v[i]*v[j] + (j == i + n) for i in 1:2*n-1 for j in i+1:2*n]
return pbw_algebra(r, strictly_upper_triangular_matrix(rel), default_ordering(r))
return pbw_algebra(r, strictly_upper_triangular_matrix(rel), default_ordering(r); check = false)
end

function weyl_algebra(
Expand Down Expand Up @@ -480,32 +488,16 @@ function base_ring(a::PBWAlgIdeal)
return a.basering
end

function _any_gens(a::PBWAlgIdeal)
return isdefined(a, :sdata) ? a.sdata : a.sopdata
end

function _set_opdata!(a::PBWAlgIdeal{D}, b::Singular.sideal) where D
@assert !isdefined(a, :sdata)
@assert !isdefined(a, :gb)
@assert !isdefined(a, :opgb)
b.isTwoSided = (D == 0)
a.sopdata = b
return a
end

function ngens(a::PBWAlgIdeal)
# this assumes a.sdata and a.sopdata have the same number of gens
return ngens(_any_gens(a))
return ngens(a.sdata)
end

function gens(a::PBWAlgIdeal{D, T, S}) where {D, T, S}
_sdata_assure!(a)
R = base_ring(a)
return PBWAlgElem{T, S}[PBWAlgElem(R, x) for x in gens(a.sdata)]
end

function gen(a::PBWAlgIdeal, i::Int)
_sdata_assure!(a)
R = base_ring(a)
return PBWAlgElem(R, a.sdata[i])
end
Expand Down Expand Up @@ -597,14 +589,6 @@ function right_ideal(R::PBWAlgRing{T, S}, g::AbstractVector) where {T, S}
return PBWAlgIdeal{1, T, S}(R, i)
end

# assure a.sdata is defined
function _sdata_assure!(a::PBWAlgIdeal)
if !isdefined(a, :sdata)
R = base_ring(a)
I.sdata = _opmap(R, a.sopdata, _opposite(R))
end
end

# assure a.sopdata is defined
function _sopdata_assure!(a::PBWAlgIdeal)
if !isdefined(a, :sopdata)
Expand All @@ -614,15 +598,14 @@ function _sopdata_assure!(a::PBWAlgIdeal)
end

function groebner_assure!(a::PBWAlgIdeal)
_sdata_assure!(a)
if !isdefined(a, :gb)
a.gb = Singular.std(a.sdata)
end
end

function opgroebner_assure!(a::PBWAlgIdeal)
_sopdata_assure!(a)
if !isdefined(a, :gb)
if !isdefined(a, :opgb)
a.opgb = Singular.std(a.sopdata)
end
end
Expand All @@ -645,7 +628,7 @@ false
```
"""
function iszero(a::PBWAlgIdeal)
return iszero(_any_gens(a))
return iszero(a.sdata)
end

function _one_check(I::Singular.sideal)
Expand Down Expand Up @@ -697,10 +680,10 @@ false
```
"""
function isone(a::PBWAlgIdeal{D}) where D
if iszero(_any_gens(a))
if iszero(a.sdata)
return false
end
if _one_check(_any_gens(a))
if _one_check(a.sdata)
return true
end
if D > 0
Expand All @@ -721,24 +704,50 @@ function Base.:+(a::PBWAlgIdeal{D, T, S}, b::PBWAlgIdeal{D, T, S}) where {D, T,
return PBWAlgIdeal{D, T, S}(base_ring(a), a.sdata + b.sdata)
end

#### To be rewritten for two-sided ideals only
####@doc Markdown.doc"""
#### *(I::PBWAlgIdeal{D, T, S}, J::PBWAlgIdeal{D, T, S}) where {D, T, S}
####
####Return the product of `I` and `J`.
####"""
####function Base.:*(a::PBWAlgIdeal{D, T, S}, b::PBWAlgIdeal{D, T, S}) where {D, T, S}
#### return PBWAlgIdeal{D, T, S}(base_ring(a), a.sdata + b.sdata)
####end

####@doc Markdown.doc"""
#### ^(I::PBWAlgIdeal{D, T, S}, m::Int) where {D, T, S}
####
####Return the `m`-th power of `I`.
####"""
####function Base.:^(a::PBWAlgIdeal{D, T, S}, b::Int) where {D, T, S}
#### return PBWAlgIdeal{D, T, S}(base_ring(a), a.sdata^b)
####end

function _as_left_ideal(a::PBWAlgIdeal{D}) where D
is_left(a) || error("cannot convert to left ideal")
if D < 0
return a.sdata
else
groebner_assure!(a)
return a.gb
end
end

function _as_right_ideal(a::PBWAlgIdeal{D}) where D
is_right(a) || error("cannot convert to right ideal")
if D > 0
return a.sdata
else
opgroebner_assure!(a)
R = base_ring(a)
return _opmap(R, a.opgb, _opposite(R))
end
end

function Base.:*(a::PBWAlgIdeal{Da, T, S}, b::PBWAlgIdeal{Db, T, S}) where {Da, Db, T, S}
@assert base_ring(a) == base_ring(b)
is_left(a) && is_right(b) || throw(NotImplementedError(:*, a, b))
# Singular.jl's cartesian product is correct for left*right
return PBWAlgIdeal{0, T, S}(base_ring(a), _as_left_ideal(a)*_as_right_ideal(b))
end

function Base.:^(a::PBWAlgIdeal{D, T, S}, b::Int) where {D, T, S}
@assert b >= 0
b == 0 && return PBWAlgIdeal{D, T, S}(base_ring(a), [one(base_ring(a))])
b == 1 && return a
if D == 0
# Note: repeated mul seems better than nested squaring
res = a
while (b -= 1) > 0
res = res*a
end
return res
else
throw(NotImplementedError(:^, a, b))
end
end

@doc Markdown.doc"""
intersect(I::PBWAlgIdeal{D, T, S}, Js::PBWAlgIdeal{D, T, S}...) where {D, T, S}
Expand All @@ -760,10 +769,8 @@ function Base.intersect(a::PBWAlgIdeal{D, T, S}, b::PBWAlgIdeal{D, T, S}...) whe
@assert R === base_ring(bi)
end
if D < 0
_sdata_assure!(a)
res = a.sdata
for bi in b
_sdata_assure!(bi)
res = Singular.intersection(res, bi.sdata)
end
return PBWAlgIdeal{D, T, S}(R, res)
Expand All @@ -774,10 +781,13 @@ function Base.intersect(a::PBWAlgIdeal{D, T, S}, b::PBWAlgIdeal{D, T, S}...) whe
_sopdata_assure!(bi)
res = Singular.intersection(res, bi.sopdata)
end
return _set_opdata!(PBWAlgIdeal{D, T, S}(R), res)
return PBWAlgIdeal{D, T, S}(R, _opmap(R, res, _opposite(R)), res)
else
# two-sided algorithm missing
throw(NotImplementedError(:intersect, a, b...))
res = _as_left_ideal(a)
for bi in b
res = Singular.intersection(res, _as_left_ideal(bi))
end
return PBWAlgIdeal{D, T, S}(R, res)
end
end

Expand Down Expand Up @@ -851,7 +861,6 @@ function Base.issubset(a::PBWAlgIdeal{D, T, S}, b::PBWAlgIdeal{D, T, S}) where {
@assert base_ring(a) === base_ring(b)
if D <= 0
# Ditto comment ideal_membership
_sdata_assure!(a)
groebner_assure!(b)
return Singular.iszero(Singular.reduce(a.sdata, b.gb))
else
Expand Down
33 changes: 29 additions & 4 deletions test/Rings/PBWAlgebra-test.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
@testset "PBWAlgebra.constructor" begin
r, (x, y, z) = QQ["x", "y", "z"]
R, (x, y, z) = pbw_algebra(r, [0 x*y x*z; 0 0 y*z + 1; 0 0 0], deglex(gens(r)))
R, (x, y, z) = pbw_algebra(r, [0 x*y x*z; 0 0 y*z + 1; 0 0 0], deglex(r))

@test elem_type(R) == PBWAlgElem{fmpq, Singular.n_Q}
@test parent_type(x) == PBWAlgRing{fmpq, Singular.n_Q}
Expand All @@ -13,11 +13,16 @@
@test R[2] == y

3*x^2 + y*z == R([3, 1], [[2, 0, 0], [0, 1, 1]])

r, (x, y, z) = QQ["x", "y", "z"]
@test_throws Exception pbw_algebra(r, [0 x*y+y x*z+z; 0 0 y*z+1; 0 0 0], lex(r))
R, (x, y, z) = pbw_algebra(r, [0 x*y+y x*z+z; 0 0 y*z+1; 0 0 0], deglex(r); check = false)
@test (z*y)*x != z*(y*x)
end

@testset "PBWAlgebra.printing" begin
r, (x, y, z) = QQ["x", "y", "z"]
R, (x, y, z) = pbw_algebra(r, [0 x*y x*z; 0 0 y*z + 1; 0 0 0], deglex(gens(r)))
R, (x, y, z) = pbw_algebra(r, [0 x*y x*z; 0 0 y*z + 1; 0 0 0], deglex(r))

@test length(string(R)) > 2
@test length(string(x + y)) > 2
Expand All @@ -26,7 +31,7 @@ end

@testset "PBWAlgebra.iteration" begin
r, (x, y, z) = QQ["x", "y", "z"]
R, (x, y, z) = pbw_algebra(r, [0 x*y x*z; 0 0 y*z + 1; 0 0 0], deglex(gens(r)))
R, (x, y, z) = pbw_algebra(r, [0 x*y x*z; 0 0 y*z + 1; 0 0 0], deglex(r))

p = -((x*z*y)^6 - 1)

Expand Down Expand Up @@ -113,7 +118,7 @@ end
@test dy*dx*x in I
@test !(x*dy*dx in I)

@test_throws NotImplementedError intersect(two_sided_ideal(R, [dy]), two_sided_ideal(R, [x]))
@test is_one(intersect(two_sided_ideal(R, [dy]), two_sided_ideal(R, [x])))

I = intersect(left_ideal([dx]), left_ideal([dy]), left_ideal([x]))
@test x^2*dx == (x*dx-1)*x
Expand All @@ -127,3 +132,23 @@ end

@test intersect(left_ideal([dx])) == left_ideal([dx])
end

@testset "PBWAlgebra.ideals.multiplication" begin
r, (x, y, z) = QQ["x", "y", "z"]
R, (x, y, z) = pbw_algebra(r, [0 x*y+y x*z+z+y; 0 0 y*z; 0 0 0], lex(r))

e1 = y*z
e2 = x*y*z

@test !is_one(two_sided_ideal([e1]))
@test !is_one(two_sided_ideal([e2]))

I = two_sided_ideal([e1*e2])

@test I == left_ideal([e1])*right_ideal([e2])
@test I != two_sided_ideal([e1])*right_ideal([e2])
@test issubset(I, two_sided_ideal([e1])*right_ideal([e2]))
@test issubset(I, two_sided_ideal([e1])*two_sided_ideal([e2]))

@test I^4 == (I^2)^2
end