search

thofma / Hecke.jl

Computational algebraic number theory
Other
224 stars 64 forks source link

Stack overflow when factoring quotients of fmpq_(m)polys #1032

Closed alexjbest closed 1 year ago

alexjbest commented 1 year ago

Not sure if I'm trying to do something silly here (real example is actually multivariate, this example seems minimal: I would expect at least a helpful error message!

julia> using Oscar

julia> R, a = PolynomialRing(QQ, ["a"])
(Multivariate Polynomial Ring in a over Rational Field, fmpq_mpoly[a])

julia> a
1-element Vector{fmpq_mpoly}:
 a

julia> factor(a[1]//ZZ(1))
ERROR: StackOverflowError:
Stacktrace:
     [1] factor(x::AbstractAlgebra.Generic.Frac{fmpq_mpoly})
       @ Hecke ~/.julia/packages/Hecke/nbhqx/src/Hecke.jl:230
     [2] factor(a::AbstractAlgebra.Generic.Frac{fmpq_mpoly}; b::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
       @ Hecke ~/.julia/packages/Hecke/nbhqx/src/Misc/Integer.jl:721
     [3] factor(a::AbstractAlgebra.Generic.Frac{fmpq_mpoly})
       @ Hecke ~/.julia/packages/Hecke/nbhqx/src/Misc/Integer.jl:721--- the last 3 lines are repeated 26660 more times ---

univariate also fails

julia> R, a = PolynomialRing(QQ)
(Univariate Polynomial Ring in x over Rational Field, x)

julia> factor(a//ZZ(1))
ERROR: StackOverflowError:
Stacktrace:
     [1] factor(a::AbstractAlgebra.Generic.Frac{fmpq_poly})
       @ Hecke ~/.julia/packages/Hecke/nbhqx/src/Misc/Integer.jl:721
     [2] factor(x::AbstractAlgebra.Generic.Frac{fmpq_poly})
       @ Hecke ~/.julia/packages/Hecke/nbhqx/src/Hecke.jl:230
     [3] factor(a::AbstractAlgebra.Generic.Frac{fmpq_poly}; b::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
       @ Hecke ~/.julia/packages/Hecke/nbhqx/src/Misc/Integer.jl:721--- the last 3 lines are repeated 26660 more times ---
thofma commented 1 year ago

Ah, we are missing the factor(x::FieldElem) fallback. Will

julia> factor(a[1]//ZZ(1))
a

make you more happy?

fieker commented 1 year ago

Apart from this being silly - what do you expect to see?

julia> a[1]//ZZ(1) a

julia> typeof(ans) AbstractAlgebra.Generic.Frac{QQMPolyRingElem}

in general, Generic.Frac does not have any factoring methods attached... Unfortunately, Hecke tries to generically dispatch and that recurses.

Was the plan/ intend to get s.th. like a rational factored with negative exponents? We can add that...

On Wed, Apr 05, 2023 at 05:12:06PM -0700, Alex J Best wrote:

Not sure if I'm trying to do something silly here (real example is actually multivariate, this example seems minimal: I would expect at least a helpful error message!

julia> using Oscar

julia> R, a = PolynomialRing(QQ, ["a"])
(Multivariate Polynomial Ring in a over Rational Field, fmpq_mpoly[a])

julia> a
1-element Vector{fmpq_mpoly}:
 a

julia> factor(a[1]//ZZ(1))
ERROR: StackOverflowError:
Stacktrace:
     [1] factor(x::AbstractAlgebra.Generic.Frac{fmpq_mpoly})
       @ Hecke ~/.julia/packages/Hecke/nbhqx/src/Hecke.jl:230
     [2] factor(a::AbstractAlgebra.Generic.Frac{fmpq_mpoly}; b::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
       @ Hecke ~/.julia/packages/Hecke/nbhqx/src/Misc/Integer.jl:721
     [3] factor(a::AbstractAlgebra.Generic.Frac{fmpq_mpoly})
       @ Hecke ~/.julia/packages/Hecke/nbhqx/src/Misc/Integer.jl:721--- the last 3 lines are repeated 26660 more times ---

univariate also fails

julia> R, a = PolynomialRing(QQ)
(Univariate Polynomial Ring in x over Rational Field, x)

julia> factor(a//ZZ(1))
ERROR: StackOverflowError:
Stacktrace:
     [1] factor(a::AbstractAlgebra.Generic.Frac{fmpq_poly})
       @ Hecke ~/.julia/packages/Hecke/nbhqx/src/Misc/Integer.jl:721
     [2] factor(x::AbstractAlgebra.Generic.Frac{fmpq_poly})
       @ Hecke ~/.julia/packages/Hecke/nbhqx/src/Hecke.jl:230
     [3] factor(a::AbstractAlgebra.Generic.Frac{fmpq_poly}; b::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
       @ Hecke ~/.julia/packages/Hecke/nbhqx/src/Misc/Integer.jl:721--- the last 3 lines are repeated 26660 more times ---

-- Reply to this email directly or view it on GitHub: https://github.com/thofma/Hecke.jl/issues/1032 You are receiving this because you are subscribed to this thread.

Message ID: @.***>

alexjbest commented 1 year ago

Well I guess I'm a bit confused that dividing a polynomial over the rationals by an integer gives me something in the fraction field in the first place, is that unavoidable for type reasons? For convenience a factoring of the numerator and denominator would be very helpful as at least the output would look nice even if the type wasn't what I was expecting.

julia> factor(a[1]//ZZ(1))
a

makes me more happy yes! Though I would love to know what your proposal does on less trivial inputs :)

thofma commented 1 year ago

The rule is that // always creates something in the (total) field of fractions. Unless you are in a field, then you stay in the field. If you want to do an (exact) division in the ring, you can do x/y:

julia> factor(a[1]^2/ZZ(1))
1 * a^2

The proposal is that

Hecke.factor(x::FieldElem) = Hecke.Fac(x, Dict{typeof(x), Int}())

thus over a field one always gets the trivial factorization consisting of the element itself (since it is a unit).

alexjbest commented 1 year ago

Yeah that makes sense. Is there a binary operation that lets me divide elements of Q algebras by integers? I guess I'm confused because I first tried these

julia> a/ ZZ(27)
ERROR: MethodError: no method matching /(::fmpq_mpoly, ::fmpz)
Closest candidates are:
  /(::ChainRulesCore.AbstractZero, ::Any) at ~/.julia/packages/ChainRulesCore/a4mIA/src/tangent_types/abstract_zero.jl:29
  /(::ChainRulesCore.AbstractThunk, ::Any) at ~/.julia/packages/ChainRulesCore/a4mIA/src/tangent_types/thunks.jl:37
  /(::arb_poly, ::fmpz) at ~/.julia/packages/Nemo/g02iz/src/arb/arb_poly.jl:324
  ...
Stacktrace:
 [1] top-level scope
   @ REPL[168]:1

julia> a / QQ(27)
ERROR: MethodError: no method matching /(::fmpq_mpoly, ::fmpq)
Closest candidates are:
  /(::ChainRulesCore.AbstractZero, ::Any) at ~/.julia/packages/ChainRulesCore/a4mIA/src/tangent_types/abstract_zero.jl:29
  /(::ChainRulesCore.AbstractThunk, ::Any) at ~/.julia/packages/ChainRulesCore/a4mIA/src/tangent_types/thunks.jl:37
  /(::arb_poly, ::fmpq) at ~/.julia/packages/Nemo/g02iz/src/arb/arb_poly.jl:324
  ...
Stacktrace:
 [1] top-level scope
   @ REPL[169]:1

julia> a / R(27)
ERROR: MethodError: no method matching /(::fmpq_mpoly, ::fmpq_mpoly)
Closest candidates are:
  /(::ChainRulesCore.AbstractZero, ::Any) at ~/.julia/packages/ChainRulesCore/a4mIA/src/tangent_types/abstract_zero.jl:29
  /(::ChainRulesCore.AbstractThunk, ::Any) at ~/.julia/packages/ChainRulesCore/a4mIA/src/tangent_types/thunks.jl:37
  /(::Any, ::ChainRulesCore.AbstractThunk) at ~/.julia/packages/ChainRulesCore/a4mIA/src/tangent_types/thunks.jl:38
  ...
Stacktrace:
 [1] top-level scope
   @ REPL[170]:1

all of which I would consider an exact ring division, but none of them work.

I understand that it is technically correct that the only factorization over a field is the trivial one, but when the field is the fraction field of a UFD it seems that factoring the numerator and denominator is often what the user would want. It would be nice if there was a function for this but I'm not sure if that would break the assumptions of library code if it were called "factor" (but perhaps a log message could be emmitted if the user calls factor on a field element telling them to try this other function) Anyway your call on whether this is worth thinking about, I can always factor(numerator(blah)) after all

thofma commented 1 year ago

If we always do what the user wants, we end up with an inconsistent mess. The function factor must do something mathematical, so it factors into irreducible elements. It will definitely break library code if we do random things for random objects. It would also mean that isirreducible does the mathematical wrong thing.

Which version do you have? What is the output of ]st -m?

thofma commented 1 year ago

Here is what you probably want:

julia> factor(QQ(9), ZZ)
1 * 3^2

julia> Qx, x = QQ["x"];

julia> factor(x//2, Qx)
(1//2) * x

julia> factor((x + 1)//(x - 2)^2, Qx)
1 * (x - 2)^-2 * (x + 1)

julia> Qt, t = QQ["t1", "t2"]
(Multivariate polynomial ring in 2 variables over QQ, QQMPolyRingElem[t1, t2])

julia> factor((t[2] + 1)//(t[1] - 2)^2, Qt)
1 * (t1 - 2)^-2 * (t2 + 1)
alexjbest commented 1 year ago

Ok thats great, I can remember to supply the ring to factor.

I would still like to be able to do exact division of fmpq_mpolys though, I think I'm on latest. Maybe this should move to nemo though?

julia> using Oscar
 -----    -----    -----      -      -----   
|     |  |     |  |     |    | |    |     |  
|     |  |        |         |   |   |     |  
|     |   -----   |        |     |  |-----   
|     |        |  |        |-----|  |   |    
|     |  |     |  |     |  |     |  |    |   
 -----    -----    -----   -     -  -     -  

...combining (and extending) ANTIC, GAP, Polymake and Singular
Version 0.11.3 ... 
 ... which comes with absolutely no warranty whatsoever
Type: '?Oscar' for more information
(c) 2019-2023 by The OSCAR Development Team

julia> R, (a, b) = PolynomialRing(QQ, ["a", "b"])
(Multivariate Polynomial Ring in a, b over Rational Field, fmpq_mpoly[a, b])

julia> a / ZZ(2)
ERROR: MethodError: no method matching /(::fmpq_mpoly, ::fmpz)
Closest candidates are:
  /(::ChainRulesCore.AbstractZero, ::Any) at ~/.julia/packages/ChainRulesCore/a4mIA/src/tangent_types/abstract_zero.jl:29
  /(::ChainRulesCore.AbstractThunk, ::Any) at ~/.julia/packages/ChainRulesCore/a4mIA/src/tangent_types/thunks.jl:37
  /(::acb_poly, ::fmpz) at ~/.julia/packages/Nemo/g02iz/src/arb/acb_poly.jl:332
  ...
Stacktrace:
 [1] top-level scope
   @ REPL[4]:1

julia> a / QQ(2)
ERROR: MethodError: no method matching /(::fmpq_mpoly, ::fmpq)
Closest candidates are:
  /(::ChainRulesCore.AbstractZero, ::Any) at ~/.julia/packages/ChainRulesCore/a4mIA/src/tangent_types/abstract_zero.jl:29
  /(::ChainRulesCore.AbstractThunk, ::Any) at ~/.julia/packages/ChainRulesCore/a4mIA/src/tangent_types/thunks.jl:37
  /(::acb_poly, ::fmpq) at ~/.julia/packages/Nemo/g02iz/src/arb/acb_poly.jl:332
  ...
Stacktrace:
 [1] top-level scope
   @ REPL[5]:1

julia> a / R(2)
ERROR: MethodError: no method matching /(::fmpq_mpoly, ::fmpq_mpoly)
Closest candidates are:
  /(::ChainRulesCore.AbstractZero, ::Any) at ~/.julia/packages/ChainRulesCore/a4mIA/src/tangent_types/abstract_zero.jl:29
  /(::ChainRulesCore.AbstractThunk, ::Any) at ~/.julia/packages/ChainRulesCore/a4mIA/src/tangent_types/thunks.jl:37
  /(::Any, ::ChainRulesCore.AbstractThunk) at ~/.julia/packages/ChainRulesCore/a4mIA/src/tangent_types/thunks.jl:38
  ...
Stacktrace:
 [1] top-level scope
   @ REPL[6]:1

julia> ]st -m
ERROR: syntax: unexpected "]"
Stacktrace:
 [1] top-level scope
   @ none:1

(@v1.8) pkg> st -m
Status `~/.julia/environments/v1.8/Manifest.toml`
⌅ [c3fe647b] AbstractAlgebra v0.27.10
  [66b61cbe] AlgebraicSolving v0.3.1
  [b99e7846] BinaryProvider v0.5.10
  [f01c122e] BinaryWrappers v0.1.2
  [d360d2e6] ChainRulesCore v1.15.7
  [9e997f8a] ChangesOfVariables v0.1.6
  [18c793ea] Coleman v0.1.0 `~/.julia/dev/Coleman`
  [34da2185] Compat v4.6.1
  [1f15a43c] CxxWrap v0.13.3
  [55939f99] DecFP v1.3.1
  [ffbed154] DocStringExtensions v0.9.3
  [c863536a] GAP v0.9.4
  [d5909c97] GroupsCore v0.4.0
⌅ [3e1990a7] Hecke v0.17.5
  [3587e190] InverseFunctions v0.1.8
  [92d709cd] IrrationalConstants v0.2.2
  [692b3bcd] JLLWrappers v1.4.1
  [682c06a0] JSON v0.21.3
⌃ [472f376f] LoadFlint v0.1.3
  [2ab3a3ac] LogExpFunctions v0.3.23
  [1914dd2f] MacroTools v0.5.10
⌅ [4fe8b98c] Mongoc v0.8.1
⌅ [2edaba10] Nemo v0.32.7
  [f1435218] Oscar v0.11.3
  [69de0a69] Parsers v2.5.8
  [d720cf60] Polymake v0.9.0
  [21216c6a] Preferences v1.3.0
  [fb686558] RandomExtensions v0.4.3
  [3cdcf5f2] RecipesBase v1.3.3
  [ae029012] Requires v1.3.0
  [6c6a2e73] Scratch v1.2.0
⌅ [bcd08a7b] Singular v0.16.1
  [66db9d55] SnoopPrecompile v1.0.3
  [276daf66] SpecialFunctions v2.2.0
  [ae81ac8f] ASL_jll v0.1.3+0
  [e21ec000] Antic_jll v0.201.500+0
  [d9960996] Arb_jll v200.2300.0+0
  [6e34b625] Bzip2_jll v1.0.8+0
  [fcfa6d1b] Calcium_jll v0.401.100+0
  [47200ebd] DecFP_jll v2.0.3+1
  [e134572f] FLINT_jll v200.900.5+0
  [5cd7a574] GAP_jll v400.1200.200+0
  [de1ad85e] GAP_lib_jll v400.1201.200+0
⌅ [ba154793] GAP_pkg_juliainterface_jll v0.800.201+2
  [e8aa6df9] GLPK_jll v5.0.1+0
  [9cc047cb] Ipopt_jll v300.1400.1000+0
  [1d63c593] LLVMOpenMP_jll v15.0.4+0
  [dd4b983a] LZO_jll v2.10.1+0
  [d00139f3] METIS_jll v5.1.2+0
  [d7ed1dd3] MUMPS_seq_jll v500.500.101+0
  [90100e71] MongoC_jll v1.19.1+0
  [68e3532b] Ncurses_jll v6.2.0+0
  [76642167] Ninja_jll v1.11.1+0
⌅ [656ef2d0] OpenBLAS32_jll v0.3.17+0
  [458c3c95] OpenSSL_jll v1.1.20+0
  [efe28fd5] OpenSpecFun_jll v0.5.5+0
  [80dd9cbb] PPL_jll v1.2.1+0
  [83958c19] Perl_jll v5.34.0+2
  [05236dd9] Readline_jll v8.1.1+1
  [e5ac4fe4] SCIP_jll v800.0.301+0
⌅ [43d676ae] Singular_jll v403.101.500+0
  [36f60fef] TOPCOM_jll v0.17.8+0
  [3161d3a3] Zstd_jll v1.5.4+0
  [508c9074] bliss_jll v0.77.0+1
⌅ [28df3c45] boost_jll v1.76.0+1
  [f07e07eb] cddlib_jll v0.94.13+0
  [5558cf25] cohomCalg_jll v0.32.0+0
  [1493ae25] lib4ti2_jll v1.6.10+0
  [3eaa8342] libcxxwrap_julia_jll v0.9.4+0
  [4d8266f6] libpolymake_julia_jll v0.9.0+0
⌅ [ae4fbd8f] libsingular_julia_jll v0.29.0+0
  [3873f7d0] lrslib_jll v0.3.3+0
  [6d01cc9a] msolve_jll v0.4.9+0
  [55c6dc9b] nauty_jll v2.6.13+0
  [6690c6e9] normaliz_jll v300.900.301+0
  [7c209550] polymake_jll v400.900.0+0
  [fe1e1685] snappy_jll v1.1.9+1
  [0dad84c5] ArgTools v1.1.1
  [56f22d72] Artifacts
  [2a0f44e3] Base64
  [ade2ca70] Dates
  [8ba89e20] Distributed
  [f43a241f] Downloads v1.6.0
  [7b1f6079] FileWatching
  [b77e0a4c] InteractiveUtils
  [4af54fe1] LazyArtifacts
  [b27032c2] LibCURL v0.6.3
  [76f85450] LibGit2
  [8f399da3] Libdl
  [37e2e46d] LinearAlgebra
  [56ddb016] Logging
  [d6f4376e] Markdown
  [a63ad114] Mmap
  [ca575930] NetworkOptions v1.2.0
  [44cfe95a] Pkg v1.8.0
  [de0858da] Printf
  [9abbd945] Profile
  [3fa0cd96] REPL
  [9a3f8284] Random
  [ea8e919c] SHA v0.7.0
  [9e88b42a] Serialization
  [6462fe0b] Sockets
  [2f01184e] SparseArrays
  [10745b16] Statistics
  [fa267f1f] TOML v1.0.0
  [a4e569a6] Tar v1.10.1
  [8dfed614] Test
  [cf7118a7] UUIDs
  [4ec0a83e] Unicode
  [e66e0078] CompilerSupportLibraries_jll v1.0.1+0
  [781609d7] GMP_jll v6.2.1+2
  [deac9b47] LibCURL_jll v7.84.0+0
  [29816b5a] LibSSH2_jll v1.10.2+0
  [3a97d323] MPFR_jll v4.1.1+1
  [c8ffd9c3] MbedTLS_jll v2.28.0+0
  [14a3606d] MozillaCACerts_jll v2022.2.1
  [4536629a] OpenBLAS_jll v0.3.20+0
  [05823500] OpenLibm_jll v0.8.1+0
  [83775a58] Zlib_jll v1.2.12+3
  [8e850b90] libblastrampoline_jll v5.1.1+0
  [8e850ede] nghttp2_jll v1.48.0+0
  [3f19e933] p7zip_jll v17.4.0+0
Info Packages marked with ⌃ and ⌅ have new versions available, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated -m`
thofma commented 1 year ago

Ah yes, the Oscar is recent, but all the other packages have been updated since then. There is no Oscar version which includes the most recent AA version yet (which includes the change to / that you are missing). At the moment you have to do divexact(a[1], 2) instead of a[1]/2.

alexjbest commented 1 year ago

Thanks! I'll leave this open to track the stackoverflow which should be the identity function. But I'm happy I can do everything I need now :)

alexjbest commented 1 year ago

I see you opened https://github.com/Nemocas/AbstractAlgebra.jl/pull/1319 so I'll close this now

thofma commented 1 year ago

Thanks for your feedback, much appreciated.