Skip to content

Commit

Permalink
feat: adjust to factor
Browse files Browse the repository at this point in the history
  • Loading branch information
@thofma
thofma committed Sep 18, 2024
1 parent 5f8c496 commit a1cfed2
Show file tree
Hide file tree
Showing 5 changed files with 29 additions and 270 deletions.
4 changes: 2 additions & 2 deletions Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand All @@ -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"
Expand Down
8 changes: 4 additions & 4 deletions src/Hecke.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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!,]

Expand All @@ -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()
Expand Down Expand Up @@ -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()
Expand Down
177 changes: 1 addition & 176 deletions src/Misc/Integer.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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"""
Expand Down
106 changes: 22 additions & 84 deletions src/Misc/coprime.jl
View file
Original file line number Diff line number Diff line change
@@ -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)
Expand Down Expand Up @@ -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)
4 changes: 0 additions & 4 deletions src/NumField/NfAbs/MPolyFactor.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down

0 comments on commit a1cfed2

Please sign in to comment.