diff --git a/Project.toml b/Project.toml index e566d70204..e090e3d3de 100644 --- a/Project.toml +++ b/Project.toml @@ -28,7 +28,7 @@ GAPExt = "GAP" PolymakeExt = "Polymake" [compat] -AbstractAlgebra = "^0.42.0" +AbstractAlgebra = "^0.43.0" Dates = "1.6" Distributed = "1.6" GAP = "0.9.6, 0.10, 0.11" @@ -37,7 +37,7 @@ LazyArtifacts = "1.6" Libdl = "1.6" LinearAlgebra = "1.6" Markdown = "1.6" -Nemo = "^0.46.0" +Nemo = "^0.47.0" Pkg = "1.6" Polymake = "0.10, 0.11" Printf = "1.6" diff --git a/src/Hecke.jl b/src/Hecke.jl index 0885e2e87a..644089a123 100644 --- a/src/Hecke.jl +++ b/src/Hecke.jl @@ -77,7 +77,7 @@ import AbstractAlgebra: get_cached!, @alias import AbstractAlgebra: pretty, Lowercase, LowercaseOff, Indent, Dedent, terse, is_terse -import AbstractAlgebra: Solve +import AbstractAlgebra: Solve, coprime_base_steel import LinearAlgebra: dot, nullspace, rank, ishermitian @@ -94,7 +94,7 @@ import Nemo import Pkg -exclude = [:Nemo, :AbstractAlgebra, :RealNumberField, :zz, :qq, :factor, :call, +exclude = [:Nemo, :AbstractAlgebra, :RealNumberField, :zz, :qq, :call, :factors, :parseint, :strongequal, :window, :xgcd, :rows, :cols, :set_entry!,] @@ -110,7 +110,7 @@ import Nemo: acb_struct, Ring, Group, Field, zzModRing, zzModRingElem, arf_struc FpField, acb_vec, array, acb_vec_clear, force_coerce, force_op, fmpz_mod_ctx_struct, divisors, is_zero_entry, IntegerUnion, remove!, valuation!, is_cyclo_type, is_embedded, is_maxreal_type, - mat_entry_ptr + mat_entry_ptr, factor_trial_range AbstractAlgebra.@include_deprecated_bindings() Nemo.@include_deprecated_bindings() @@ -164,7 +164,7 @@ function __init__() global R = _RealRing() - global flint_rand_ctx = flint_rand_state() + #global flint_rand_ctx = flint_rand_state() resize!(_RealRings, Threads.nthreads()) for i in 1:Threads.nthreads() diff --git a/src/Misc/Integer.jl b/src/Misc/Integer.jl index f1d4b5fc37..2b760599e7 100644 --- a/src/Misc/Integer.jl +++ b/src/Misc/Integer.jl @@ -172,181 +172,7 @@ function is_prime_power(n::Integer) return is_prime_power(ZZRingElem(n)) end -################################################################################ -# random and factor -################################################################################ - -factor(a...; b...) = Nemo.factor(a...; b...) - -factor(a::Integer) = factor(ZZRingElem(a)) - -mutable struct flint_rand_ctx_t - a::Ptr{Nothing} - function flint_rand_ctx_t() - return new() - end -end - -function show(io::IO, A::flint_rand_ctx_t) - println(io, "Flint random state") -end - -function flint_rand_state() - A = flint_rand_ctx_t() - A.a = ccall((:flint_rand_alloc, libflint), Ptr{Nothing}, (Int,), 1) - ccall((:flint_randinit, libflint), Nothing, (Ptr{Nothing},), A.a) - - function clean_rand_state(A::flint_rand_ctx_t) - ccall((:flint_randclear, libflint), Nothing, (Ptr{Nothing},), A.a) - ccall((:flint_rand_free, libflint), Nothing, (Ptr{Nothing},), A.a) - nothing - end - finalizer(clean_rand_state, A) - return A -end - -global flint_rand_ctx - -function ecm(a::ZZRingElem, B1::UInt, B2::UInt, ncrv::UInt, rnd=flint_rand_ctx) - f = ZZRingElem() - r = ccall((:fmpz_factor_ecm, libflint), Int32, (Ref{ZZRingElem}, UInt, UInt, UInt, Ptr{Nothing}, Ref{ZZRingElem}), f, ncrv, B1, B2, rnd.a, a) - return r, f -end - -function ecm(a::ZZRingElem, B1::Int, B2::Int, ncrv::Int, rnd=flint_rand_ctx) - return ecm(a, UInt(B1), UInt(B2), UInt(ncrv), rnd) -end - -#data from http://www.mersennewiki.org/index.php/Elliptic_Curve_Method -const B1 = [2, 11, 50, 250, 1000, 3000, 11000, 43000, 110000, 260000, 850000, 2900000]; -const nC = [25, 90, 300, 700, 1800, 5100, 10600, 19300, 49000, 124000, 210000, 340000]; - -function ecm(a::ZZRingElem, max_digits::Int=div(ndigits(a), 3), rnd=flint_rand_ctx) - n = ndigits(a, 10) - B1s = 15 - - i = 1 - s = max(div(max_digits - 10, 5), 1) - #i = s = max(i, s) - while i <= s - e, f = ecm(a, B1[i] * 1000, B1[i] * 1000 * 100, nC[i], rnd) - if e != 0 - return (e, f) - end - i += 1 - if i > length(B1) - return (e, f) - end - end - return (Int32(0), a) -end - -function factor_trial_range(N::ZZRingElem, start::Int=0, np::Int=10^5) - F = Nemo.fmpz_factor() - ccall((:fmpz_factor_trial_range, libflint), Nothing, (Ref{Nemo.fmpz_factor}, Ref{ZZRingElem}, UInt, UInt), F, N, start, np) - res = Dict{ZZRingElem,Int}() - for i in 1:F.num - z = ZZRingElem() - ccall((:fmpz_factor_get_fmpz, libflint), Nothing, - (Ref{ZZRingElem}, Ref{Nemo.fmpz_factor}, Int), z, F, i - 1) - res[z] = unsafe_load(F.exp, i) - end - return res, canonical_unit(N) -end - -const big_primes = ZZRingElem[] - -function factor(N::ZZRingElem) - if iszero(N) - throw(ArgumentError("Argument is not non-zero")) - end - N_in = N - global big_primes - r, c = factor_trial_range(N) - for (p, v) = r - N = divexact(N, p^v) - end - if is_unit(N) - @assert N == c - return Nemo.Fac(c, r) - end - N *= c - @assert N > 0 - - for p = big_primes - v, N = remove(N, p) - if v > 0 - @assert !haskey(r, p) - r[p] = v - end - end - factor_insert!(r, N) - for p = keys(r) - if nbits(p) > 60 && !(p in big_primes) - push!(big_primes, p) - end - end - return Nemo.Fac(c, r) -end - -function factor_insert!(r::Dict{ZZRingElem,Int}, N::ZZRingElem, scale::Int=1) - #assumes N to be positive - # no small divisors - # no big_primes - if isone(N) - return r - end - fac, N = is_perfect_power_with_data(N) - if fac > 1 - return factor_insert!(r, N, fac * scale) - end - if is_prime(N) - @assert !haskey(r, N) - r[N] = scale - return r - end - if ndigits(N) < 60 - s = Nemo.factor(N) #MPQS - for (p, k) in s - if haskey(r, p) - r[p] += k * scale - else - r[p] = k * scale - end - end - return r - end - - e, f = ecm(N) - if e == 0 - s = Nemo.factor(N) - for (p, k) in s - if haskey(r, p) - r[p] += k * scale - else - r[p] = k * scale - end - end - return r - end - cp = coprime_base([N, f]) - for i = cp - factor_insert!(r, i, scale * valuation(N, i)) - end - return r -end - -#TODO: problem(s) -# Nemo.factor = mpqs is hopeless if > n digits, but asymptotically and practically -# faster than ecm. -# ecm is much better if there are "small" factors. -# p-1 and p+1 methods are missing -# so probably -# if n is small enough -> Nemo -# if n is too large: ecm -# otherwise -# need ecm to find small factors -# then recurse... +###################################################################################### function _factors_trial_division(n::ZZRingElem, np::Int=10^5) res, u = factor_trial_range(n, 0, np) @@ -356,7 +182,6 @@ function _factors_trial_division(n::ZZRingElem, np::Int=10^5) n = divexact(n, p^v) end return factors, n - end @doc raw""" diff --git a/src/Misc/coprime.jl b/src/Misc/coprime.jl index f70e854c8b..b632f18c10 100644 --- a/src/Misc/coprime.jl +++ b/src/Misc/coprime.jl @@ -1,3 +1,25 @@ + +##implemented +# Bernstein: asymptotically fastest, linear in the total input size +# pointless for small ints as it requires intermediate numbers to be +# huge +# Bach/Shallit/???: not too bad, source Epure's Masters thesis +# can operate on Int types as no intermediate is larger than the input +# Steel: a misnomer: similar to Magma, basically implements a memory +# optimised version of Bach +# faster than Magma on +# > I := [Random(1, 10000) * Random(1, 10000) : x in [1..10000]]; +# > time c := CoprimeBasis(I); +# julia> I = [ZZRingElem(rand(1:10000))*rand(1:10000) for i in 1:10000]; +# +# experimentally, unless the input is enormous, Steel wins +# on smallish input Bach is better than Bernstein, on larger this +# changes +# +# needs +# isone, gcd_into!, divexact!, copy +# (some more for Bernstein: FactorBase, gcd, divexact) + import Nemo.isone, Nemo.divexact, Base.copy function gcd_into!(a::ZZRingElem, b::ZZRingElem, c::ZZRingElem) @@ -423,87 +445,3 @@ function coprime_base_bernstein(S::Vector{E}) where E return merge_bernstein(P1, P2) end -function augment_steel(S::Vector{E}, a::E, start::Int = 1) where E - i = start - if is_unit(a) - return S - end - - g = a - new = true - - while i<=length(S) && !isone(a) - if new - g = gcd(S[i], a) - new = false - else - g = gcd_into!(g, S[i], a) - end - if is_unit(g) - i += 1 - continue - end - si = divexact(S[i], g) - a = divexact(a, g) - if is_unit(si) # g = S[i] and S[i] | a - continue - end - S[i] = si - if is_unit(a) # g = a and a | S[i] - a = copy(g) - continue - end - augment_steel(S, copy(g), i) - continue - end - if !is_unit(a) - push!(S, a) - end - - return S; -end - -function coprime_base_steel(S::Vector{E}) where E - @assert !isempty(S) - T = Array{E}(undef, 1) - T[1] = S[1] - for i=2:length(S) - augment_steel(T, S[i]) - end - return T -end - -##implemented -# Bernstein: asymptotically fastest, linear in the total input size -# pointless for small ints as it requires intermediate numbers to be -# huge -# Bach/Shallit/???: not too bad, source Epure's Masters thesis -# can operate on Int types as no intermediate is larger than the input -# Steel: a misnomer: similar to Magma, basically implements a memory -# optimised version of Bach -# faster than Magma on -# > I := [Random(1, 10000) * Random(1, 10000) : x in [1..10000]]; -# > time c := CoprimeBasis(I); -# julia> I = [ZZRingElem(rand(1:10000))*rand(1:10000) for i in 1:10000]; -# -# experimentally, unless the input is enormous, Steel wins -# on smallish input Bach is better than Bernstein, on larger this -# changes -# -# needs -# isone, gcd_into!, divexact!, copy -# (some more for Bernstein: FactorBase, gcd, divexact) - -@doc raw""" - coprime_base{E}(S::Vector{E}) -> Vector{E} - -Returns a coprime base for $S$, i.e. the resulting array contains pairwise coprime objects that multiplicatively generate the same set as the input array. -""" -coprime_base(x) = coprime_base_steel(x) - -@doc raw""" - coprime_base_insert{E}(S::Vector{E}, a::E) -> Vector{E} - -Given a coprime array $S$, insert a new element, i.e. find a coprime base for `push(S, a)`. -""" -coprime_base_insert(S, a) = augment_steel(S, a) diff --git a/src/NumField/NfAbs/MPolyFactor.jl b/src/NumField/NfAbs/MPolyFactor.jl index 0d830e9ed1..948b965731 100644 --- a/src/NumField/NfAbs/MPolyFactor.jl +++ b/src/NumField/NfAbs/MPolyFactor.jl @@ -25,10 +25,6 @@ function AbstractAlgebra.factor_squarefree(a::Generic.MPoly{<:NumFieldElem}) return _doit(a, factor_squarefree) end -function AbstractAlgebra.factor(a::Generic.Poly{AbsSimpleNumFieldElem}) - return Hecke.factor(a) -end - function AbstractAlgebra.MPolyFactor.mfactor_choose_eval_points!( alphas::Vector, A::E,