diff --git a/src/EllCrv/EllCrv.jl b/src/EllCrv/EllCrv.jl index 7ffc73dd08..002f2cd8b9 100644 --- a/src/EllCrv/EllCrv.jl +++ b/src/EllCrv/EllCrv.jl @@ -164,11 +164,11 @@ end function Base.getindex(P::EllCrvPt, i::Int) @req 1 <= i <= 3 "Index must be 1, 2 or 3" - + K = base_field(parent(P)) - + if is_infinite(P) - if i == 1 + if i == 1 return zero(K) elseif i == 2 return one(K) @@ -236,14 +236,14 @@ function elliptic_curve(x::Vector{Rational{T}}; check::Bool = true) where {T <: return elliptic_curve(QQFieldElem[QQ(z) for z in x], check = check) end -# A constructor to create an elliptic curve from a bivariate polynomial. -# One can specify how to interpret the polynomial via the second and the +# A constructor to create an elliptic curve from a bivariate polynomial. +# One can specify how to interpret the polynomial via the second and the # third argument. @doc raw""" elliptic_curve(f::MPolyRingElem, x::MPolyRingElem, y::MPolyRingElem) -> EllCrv -Construct an elliptic curve from a bivariate polynomial `f` in long Weierstrass form. -The second and third argument specify variables of the `parent` of `f` so that +Construct an elliptic curve from a bivariate polynomial `f` in long Weierstrass form. +The second and third argument specify variables of the `parent` of `f` so that ``f = c ⋅ (-y² + x³ - a₁ ⋅ xy + a₂ ⋅ x² - a₃ ⋅ y + a₄ ⋅ x + a₆)``. """ function elliptic_curve(f::MPolyRingElem, x::MPolyRingElem, y::MPolyRingElem) @@ -624,10 +624,10 @@ y^2 = x^3 + x + 2 ``` """ function (E::EllCrv{T})(coords::Vector{S}; check::Bool = true) where {S, T} - if !(2 <= length(coords) <= 3) + if !(2 <= length(coords) <= 3) error("Points need to be given in either affine coordinates (x, y) or projective coordinates (x, y, z)") end - + if length(coords) == 3 if coords[1] == 0 && coords[3] == 0 if coords[2] != 0 @@ -636,7 +636,7 @@ function (E::EllCrv{T})(coords::Vector{S}; check::Bool = true) where {S, T} error("The triple [0: 0: 0] does not define a point in projective space.") end end - coords = [coords[1]//coords[3], coords[2]//coords[3]] + coords = [coords[1]//coords[3], coords[2]//coords[3]] end if S === T parent(coords[1]) != base_field(E) && @@ -1122,18 +1122,18 @@ function multiplication_by_m_y_coord(E::EllCrv, m::S, x = polynomial_ring(base_f psi_m_univ = division_polynomial_univariate(E, m, x)[2] psi_mplus = division_polynomial(E, m+1, x, y) psi_mmin2 = division_polynomial(E, m-2, x, y) - - if p == 3 && j_invariant(E) != 0 - if iseven(m) + + if p == 3 && j_invariant(E) != 0 + if iseven(m) num = (psi_mplus2*psi_mmin^2 - psi_mmin2*psi_mplus^2) - a1*(x* psi_m_univ^2*B6 -psi_mmin*psi_mplus)*psi_m_univ denom = 2*B6^2*psi_m_univ^3 else num = (psi_mplus2*psi_mmin^2 - psi_mmin2*psi_mplus^2) - a1*(x* psi_m_univ^2 -psi_mmin*psi_mplus*B6)*psi_m_univ - denom = 4*y*psi_m_univ^3 + denom = 4*y*psi_m_univ^3 end return num//denom end - + num = (psi_mplus2*psi_mmin^2 - psi_mmin2*psi_mplus^2) diff --git a/src/EllCrv/Finite.jl b/src/EllCrv/Finite.jl index b768847a70..cc586fed4e 100644 --- a/src/EllCrv/Finite.jl +++ b/src/EllCrv/Finite.jl @@ -35,7 +35,7 @@ ################################################################################ export hasse_interval, order, order_via_exhaustive_search, order_via_bsgs, order_via_legendre, - order_via_schoof, trace_of_frobenius, rand, elem_order_bsgs, is_supersingular, + order_via_schoof, trace_of_frobenius, rand, elem_order_bsgs, is_supersingular, is_ordinary, is_probable_supersingular, supersingular_polynomial ################################################################################ @@ -924,34 +924,34 @@ Return true when the elliptic curve is supersingular. The result is proven to be """ function is_supersingular(E::EllCrv{T}) where T <: FinFieldElem K = base_field(E) - + p = characteristic(K) j = j_invariant(E) - + if j^(p^2) != j return false end - + if p<= 3 return j == 0 end - + L = Native.GF(p, 2) Lx, X = polynomial_ring(L, "X") Lxy, Y = polynomial_ring(Lx, "Y") Phi2 = X^3 + Y^3 - X^2*Y^2 + 1488*(X^2*Y + Y^2*X) - 162000*(X^2 + Y^2) + 40773375*X*Y + 8748000000*(X + Y) - 157464000000000 - + jL = _embed_into_p2(j, L) - + js = roots(Phi2(jL)) - + if length(js) < 3 return false end - + newjs = [jL, jL, jL] f = elem_type(Lx)[zero(Lx), zero(Lx), zero(Lx)] - + m = nbits(p) - 1 for k in (1 : m) for i in (1 : 3) @@ -1016,7 +1016,7 @@ function is_probable_supersingular(E::EllCrv{T}) where T <: FinFieldElem j = j_invariant(E) K = base_field(E) p = characteristic(K) - + local degj::Int if degree(K) == 1 @@ -1024,7 +1024,7 @@ function is_probable_supersingular(E::EllCrv{T}) where T <: FinFieldElem else degj = degree(minpoly(j)) end - + if degj == 1 return monte_carlo_test(E, p+1) elseif degj == 2 @@ -1036,14 +1036,14 @@ end function monte_carlo_test(E, n) E_O = infinity(E) - + for i in (1:10) P = rand(E) if n*P != E_O return false end end - + return true end @@ -1060,7 +1060,7 @@ function supersingular_polynomial(p::IntegerUnion) if p < 3 return J end - + m = divexact((p-1), 2) KXT, (X, T) = polynomial_ring(K, ["X", "T"]) H = sum([binomial(m, i)^2 *T^i for i in (0:m)]) @@ -1144,7 +1144,7 @@ Return a list of generators of the group of rational points on $E$. # Examples -```jldoctest; filter = r"Point.*" +```jldoctest; filter = r"Point.*" julia> E = elliptic_curve(GF(101, 2), [1, 2]); julia> gens(E) diff --git a/src/EllCrv/FormalGroupLaw.jl b/src/EllCrv/FormalGroupLaw.jl index 3c1bfc1d71..cf8acbf8b1 100644 --- a/src/EllCrv/FormalGroupLaw.jl +++ b/src/EllCrv/FormalGroupLaw.jl @@ -19,17 +19,17 @@ function formal_w(E::EllCrv, prec::Int = 20) k = base_field(E) kz, z = laurent_series_ring(k, prec, "z") kzw, w = polynomial_ring(kz, "w") - + a1, a2, a3, a4, a6 = a_invars(E) - + f = z^3 + a1*z*w + a2*z^2*w + a3*w^2 + a4*z*w^2 + a6*w^3 result = z^3 + a1*z*w + a2*z^2*w + a3*w^2 + a4*z*w^2 + a6*w^3 - + for i in (1:prec) result = truncate(f(result), div(10,3)) end return result(0) - + end @doc raw""" @@ -69,7 +69,7 @@ function formal_differential_form(E::EllCrv, prec::Int = 20) x = formal_x(E, prec + 1) y = formal_y(E, prec + 1) dx = derivative(x) - + return truncate(dx//(2*y + a1*x + a3), prec) end @@ -87,29 +87,29 @@ end k = base_field(E) a1, a2, a3, a4, a6 = a_invars(E) ktt, (z1, z2) = power_series_ring(k, prec + 1, ["z1", "z2"]) - + #We now compute the slope lambda of the addition in the formal group law #Following Silverman #A_n-3 = w[n] for n in 3 to infinity - #(z1^n - z2^n)//(z1 - z2) = z1^(n-1) + z1^(n-2)*z2 +.... +z2^(n-2) - #lambda = sum A_n-3*(z1^n - z2^n)//(z1 - z2) - + #(z1^n - z2^n)//(z1 - z2) = z1^(n-1) + z1^(n-2)*z2 +.... +z2^(n-2) + #lambda = sum A_n-3*(z1^n - z2^n)//(z1 - z2) + w = formal_w(E, prec + 1) - + lambda = sum([coeff(w, n)*sum([z1^m * z2^(n-m-1) for m in (0:n-1)]) for n in (3:prec +1)]) - + #Needs to simply be evaluating wz1 = sum([coeff(w, i)*z1^i for i in (0:prec +1)]) - + nu = wz1 - lambda*z1 - + z3 = -z1 - z2 - divexact(a1*lambda + a3*lambda^2 + a2*nu + 2*a4*lambda*nu + 3*a6*lambda^2*nu, 1 + a2*lambda + a4*lambda^2 + a6*lambda^3) - + inv = formal_inverse(E, prec) - + #Needs to simply be evaluating result = sum([coeff(inv, i)*z3^i for i in (0:prec +1)]) - + #Sage and Magma truncate to O(z1, z2)^20 instead of O(z1)^20 + O(z2)^20 like we do return result end @@ -138,11 +138,11 @@ function formal_isogeny(phi::Isogeny, prec::Int = 20) E = domain(phi) x = formal_x(E, prec) y = formal_y(E, prec) - + Rz = parent(x) - + maps = rational_maps(phi) - + fnum = change_base_ring(Rz, numerator(maps[1])) fdenom = change_base_ring(Rz, denominator(maps[1])) gnum = change_base_ring(Rz, numerator(maps[2])) diff --git a/src/EllCrv/Heights.jl b/src/EllCrv/Heights.jl index e285a7182f..ae57278a68 100644 --- a/src/EllCrv/Heights.jl +++ b/src/EllCrv/Heights.jl @@ -34,7 +34,7 @@ # ################################################################################ -export local_height, canonical_height, naive_height, height_pairing, +export local_height, canonical_height, naive_height, height_pairing, regulator, neron_tate_height, CPS_dvev_real, CPS_dvev_complex, CPS_non_archimedean, CPS_height_bounds, derivative, refine_alpha_bound @@ -74,35 +74,35 @@ end function naive_height_coordinate(x::nf_elem, prec::Int = 100) attempt = 1 - + K = parent(x) OK = ring_of_integers(K) - + q = K(denominator(x)) - + N = norm(ideal(OK, x) + 1*OK) - + deg = degree(K) while true R = ArbField(attempt*prec, cached = false) - + #Non-archimedean contribution result = -log(N) - + #Archimedean contribution (Mahler measure) for v in real_places(K) s = abs(evaluate(x, _embedding(v), attempt*prec)) result = result + log(max(s, one(R))) end - + for v in complex_places(K) s = abs(evaluate(x, _embedding(v), attempt*prec)) result = result + 2*log(max(s, one(R))) end - + result = result//deg - + if radiuslttwopower(result, -prec) expand!(result, -prec) @assert radiuslttwopower(result, -prec) @@ -159,7 +159,7 @@ function local_height(P::EllCrvPt{QQFieldElem}, p, prec::Int = 100) P = phi(P) p = FlintZZ(p) - + x = P[1] y = P[2] @@ -230,12 +230,12 @@ function local_height(P::EllCrvPt{nf_elem}, pIdeal::NfOrdIdl, prec::Int = 100) K = base_field(E) OK = ring_of_integers(K) F, phi = minimal_model(E, pIdeal) - + res_degree = norm(pIdeal) p = minimum(pIdeal) P = phi(P) - + x = P[1] y = P[2] @@ -350,7 +350,7 @@ function _real_height(P::EllCrvPt{QQFieldElem}, prec = 100) ) while true - R = ArbField(attempt*wprec) + R = ArbField(attempt*wprec) x = R(P[1]) y = R(P[2]) @@ -449,14 +449,14 @@ function archimedean_height(P::EllCrvPt{nf_elem}, _v::InfPlc, prec = 100) ) while true - newprec = attempt*wprec + newprec = attempt*wprec b2, b4, b6, b8 = map(t -> evaluate(t, v, newprec), get_b_integral(F)) - + _b2 = b2-12 _b4 = b4-b2+6 _b6 = b6-2*b4+b2-4 _b8 = b8-3*b6+3*b4-b2+3 - + x = evaluate(P[1], v, newprec) y = evaluate(P[2], v, newprec) @@ -518,7 +518,7 @@ end ################################################################################ @doc raw""" - neron_tate_height(P::EllCrvPt{T}, prec::Int) -> arb + neron_tate_height(P::EllCrvPt{T}, prec::Int) -> arb where T<:Union{QQFieldElem, nf_elem} Compute the Néron-Tate height (or canonical height) of a point $P$ on an @@ -574,30 +574,30 @@ function canonical_height(P::EllCrvPt{nf_elem}, prec = 100) R = ArbField(attempt*prec, cached = false) E = P.parent disc = discriminant(E) - - #d should be the norm of J where I/J = P[1]*OK is the unique decomposition + + #d should be the norm of J where I/J = P[1]*OK is the unique decomposition #of prime integer ideals d = (denominator(P[1]*OK)) h = log(d) - + for v in real_places(K) h = h + local_height(P, v, attempt*prec) end - + for v in complex_places(K) h = h + 2*local_height(P, v, attempt*prec) end - - + + plist = bad_primes(E) #Removed the divides check for p in plist h = h + local_height(P,p, attempt*prec) end - + h = h//degree(K) - + if radiuslttwopower(h, -prec) expand!(h, -prec) @assert radiuslttwopower(h, -prec) @@ -609,11 +609,11 @@ function canonical_height(P::EllCrvPt{nf_elem}, prec = 100) end @doc raw""" - height_pairing(P::EllCrvPt{T},Q::EllCrvPt{T}, prec::Int) + height_pairing(P::EllCrvPt{T},Q::EllCrvPt{T}, prec::Int) -> ArbField where T<:Union{QQFieldElem, nf_elem} Compute the height pairing of two points $P$ and $Q$ of an -elliptic curve defined over a number field. It is defined by +elliptic curve defined over a number field. It is defined by $h(P,Q) = (h(P + Q) - h(P) -h(Q))/2$ where $h$ is the canonical height. """ function height_pairing(P::EllCrvPt{T}, Q::EllCrvPt{T}, prec::Int = 100) where T<:Union{QQFieldElem, nf_elem} @@ -692,24 +692,24 @@ function CPS_height_bounds(E::EllCrv{T}) where T<:Union{QQFieldElem, nf_elem} P = bad_primes(E) K = base_field(E) d = degree(K) - + Rv = real_places(K) Cv = complex_places(K) - + dv_arch = ev_arch = zero(ArbField(prec, cached = false)) - + for v in Rv - dv, ev = CPS_dvev_real(E, v, prec) + dv, ev = CPS_dvev_real(E, v, prec) dv_arch += log(dv) ev_arch += log(ev) end - + for v in Cv - dv, ev = CPS_dvev_complex(E, v, prec) + dv, ev = CPS_dvev_complex(E, v, prec) dv_arch += 2*log(dv) ev_arch += 2*log(ev) end - + non_arch_contribution = sum([CPS_non_archimedean(E, v, prec) for v in P])//d return 1//(3*d) * dv_arch, 1//(3*d) * ev_arch + non_arch_contribution end @@ -718,13 +718,13 @@ function CPS_non_archimedean(E::EllCrv{T}, v, prec::Int = 100) where T OK = ring_of_integers(base_field(E)) Ep, K, f, c = tates_algorithm_local(E, v) k = K.ksymbol - + Rc = ArbField(prec, cached = false) - + # See Table 1 in Cremona, Prickett, Siksek Height Difference Bounds For Elliptic Curves - # over Number Fields for the values of av depending on + # over Number Fields for the values of av depending on # the Kodaira symbol and the Tamagawa number - if c == 1 + if c == 1 av = 0 elseif k > 4 m = k - 4 @@ -751,7 +751,7 @@ function CPS_non_archimedean(E::EllCrv{T}, v, prec::Int = 100) where T elseif k == -3 av = 3//2 end - + D = discriminant(E) Dv = discriminant(Ep) qv = norm(v) @@ -764,39 +764,39 @@ function CPS_dvev_real(E::EllCrv{T}, v::V, prec::Int = 100) where T where V<:Uni C = AcbField(prec) K = base_field(E) Kx, x = polynomial_ring(K, "x") - + b2, b4, b6, b8 = b_invars(E) - + f = 4*x^3 + b2*x^2 + 2*b4*x + b6 df = 12*x^2 +2*b2*x + 2*b4 g = x^4 - b4*x^2 -2*b6*x - b8 dg = 4*x^3 - 2*b4*x -2*b6 - + F = b6*x^4 + 2*b4*x^3 + b2*x^2 + 4*x G = -b8*x^4 - 2*b6*x^3 - b4*x^2 + 1 dF = 4*b6*x^3 + 6*b4*x^2 + 2*b2*x + 4 dG = -4*b8*x^3 - 6*b6*x^2 - 2*b4*x - + S = vcat(_roots(f, v, prec = prec)[2], _roots(g, v, prec = prec)[2], _roots(f + g, v, prec = prec)[2], _roots(f - g, v, prec = prec)[2], _roots(df, v, prec = prec)[2], _roots(dg, v, prec = prec)[2], Rc(1), Rc(-1)) S2 = vcat(_roots(F, v, prec = prec)[2], _roots(G, v, prec = prec)[2], _roots(F + G, v, prec = prec)[2], _roots(F - G, v, prec = prec)[2], _roots(dF, v, prec = prec)[2], _roots(dG, v, prec = prec)[2], Rc(1), Rc(-1)) - + Rx, x = polynomial_ring(Rc, "x") - + b2R, b4R, b6R, b8R = map(real, map(t -> evaluate(t, _embedding(v), prec), b_invars(E))) - + fR = 4*x^3 + b2R*x^2 + 2*b4R*x + b6R gR = x^4 - b4R*x^2 -2*b6R*x - b8R FR = b6R*x^4 + 2*b4R*x^3 + b2R*x^2 + 4*x GR = -b8R*x^4 - 2*b6R*x^3 - b4R*x^2 + 1 - + test_fg = function(x::arb) - + fx = evaluate(fR, x) - + return abs(x)<= 1 && (fx > Rc(0) || contains(fx, zero(Rc))) end - + filter!(test_fg, S) fglist = map(s -> max(abs(evaluate(fR,s)), abs(evaluate(gR,s))), S) @@ -804,49 +804,49 @@ function CPS_dvev_real(E::EllCrv{T}, v::V, prec::Int = 100) where T where V<:Uni Fx = evaluate(FR, x) return abs(x)<= 1 && (Fx > 0 ||contains(Fx, zero(Rc))) end - + filter!(test_FG, S2) - + FGlist = map(s -> max(abs(evaluate(FR,s)), abs(evaluate(GR,s))), S2) - + e_v = inv(minimum(vcat(fglist, FGlist))) d_v = inv(maximum(vcat(fglist, FGlist))) - + return d_v, e_v end function CPS_dvev_complex(E::EllCrv{T}, v::V, prec::Int = 100) where T where V<:Union{InfPlc, PosInf} - + Rc = ArbField(prec) C = AcbField(prec) K = base_field(E) Rx, x = polynomial_ring(C, "x") - + b2, b4, b6, b8 = map(t -> evaluate(t, _embedding(v), prec), b_invars(E)) - + f = 4*x^3 + b2*x^2 + 2*b4*x + b6 g = x^4 - b4*x^2 -2*b6*x - b8 F = b6*x^4 + 2*b4*x^3 + b2*x^2 + 4*x G = -b8*x^4 - 2*b6*x^3 - b4*x^2 + 1 - + E_fg = function (u::acb, eta::arb) fsum = sum([eta^i//factorial(i)*abs(derivative(f, i)(u)) for i in (1:3)]) gsum = sum([eta^i//factorial(i)*abs(derivative(g, i)(u)) for i in (1:4)]) return max(fsum, gsum) end - + E_FG = function (u::acb, eta::arb) Fsum = sum([eta^i//factorial(i)*abs(derivative(F, i)(u)) for i in (1:4)]) Gsum = sum([eta^i//factorial(i)*abs(derivative(G, i)(u)) for i in (1:4)]) return max(Fsum, Gsum) end - - + + M = [(m, n) for m in (-10:10) for n in (-10:10)] filter!(t -> t[1]^2 + t[2]^2 <= 100, M) M = map(t -> divexact((t[1]) + t[2]*onei(C), 10), M) - + H_fg = [max(abs(f(z)), abs(g(z))) for z in M] H_FG = [max(abs(F(z)), abs(G(z))) for z in M] @@ -866,16 +866,16 @@ function CPS_dvev_complex(E::EllCrv{T}, v::V, prec::Int = 100) where T where V<: end #Determines precision. Choosing a smaller mu makes the function a lot slower as - #alpha_bound converges very slowly. When I tested different values it took 0.5s for + #alpha_bound converges very slowly. When I tested different values it took 0.5s for #0.001 and 45s for 0.00001 for the same curve. Smaller mu only marginally improves the bound - #so it is probably fine. We seem to get the same results as Magma. + #so it is probably fine. We seem to get the same results as Magma. mu = Rc(0.001) a = -one(Rc) r = 2*one(Rc) - - approx_ev = inv(min(refine_alpha_bound(f, g, E_fg, mu, a, a, r, alpha_start_fg, prec), + + approx_ev = inv(min(refine_alpha_bound(f, g, E_fg, mu, a, a, r, alpha_start_fg, prec), refine_alpha_bound(F, G, E_FG, mu, a, a, r, alpha_start_FG, prec))) - approx_dv = inv(min(refine_beta_bound(f, g, E_fg, mu, a, a, r, beta_start_fg,prec), + approx_dv = inv(min(refine_beta_bound(f, g, E_fg, mu, a, a, r, beta_start_fg,prec), refine_beta_bound(F, G, E_FG, mu, a, a, r, beta_start_FG, prec))) return approx_dv, approx_ev end @@ -898,7 +898,7 @@ function refine_alpha_bound(P::PolyElem, Q::PolyElem, E, mu::arb, a::arb, b::ar if isempty(corners) && r!= 2 return alpha_bound end - + if abs(a + r//2 + (b + r//2)*i) <= 1 u = a + r//2 + (b + r//2)*i eta = r//sqrt(Rc(2)) @@ -906,20 +906,20 @@ function refine_alpha_bound(P::PolyElem, Q::PolyElem, E, mu::arb, a::arb, b::ar u = corners[1] eta = r*sqrt(Rc(2)) end - + hu = max(abs(P(u)), abs(Q(u))) if hu - E(u, eta) > alpha_bound*exp(-mu) return alpha_bound end - + alpha_bound = min(alpha_bound, hu) - + #Subdiving the initial square into four squares and computing alpha_bound there alpha_bound = refine_alpha_bound(P, Q, E, mu, a, b, r//2, alpha_bound, prec) alpha_bound = refine_alpha_bound(P, Q, E, mu, a, b + r//2, r//2, alpha_bound, prec) alpha_bound = refine_alpha_bound(P, Q, E, mu, a + r//2, b, r//2, alpha_bound, prec) alpha_bound = refine_alpha_bound(P, Q, E, mu, a + r//2, b + r//2, r//2, alpha_bound, prec) - + return alpha_bound end @@ -934,7 +934,7 @@ function refine_beta_bound(P::PolyElem, Q::PolyElem, E, mu::arb, a::arb, b::arb #We consider a square [a, a + r] x [b*i, (b + r) * i] in C corners = [a + b*i, a + r + b*i, a + (b + r)*i, a + r + (b + r)*i] filter!(t -> abs(t) < 1, corners) - + #The supremum of h(u) is attained on the boundary, so we #stop if the square doesn't intersect with the boundary. #This happens exactly when none of the corners lie either inside @@ -943,7 +943,7 @@ function refine_beta_bound(P::PolyElem, Q::PolyElem, E, mu::arb, a::arb, b::arb if (isempty(corners) || length(corners) == 4) && r!= 2 return beta_bound end - + if abs(a + r//2 + (b + r//2)*i) <= 1 u = a + r//2 + (b + r//2)*i eta = r//sqrt(Rc(2)) @@ -951,20 +951,20 @@ function refine_beta_bound(P::PolyElem, Q::PolyElem, E, mu::arb, a::arb, b::arb u = corners[1] eta = r*sqrt(Rc(2)) end - + hu = max(abs(P(u)), abs(Q(u))) - + if hu - E(u, eta) < beta_bound*exp(mu) return beta_bound end - + beta_bound = max(beta_bound, hu) - + #Subdiving the initial square into four squares and computing alpha_bound there beta_bound = refine_beta_bound(P, Q, E, mu, a, b, r//2, beta_bound, prec) beta_bound = refine_beta_bound(P, Q, E, mu, a, b + r//2, r//2, beta_bound, prec) beta_bound = refine_beta_bound(P, Q, E, mu, a + r//2, b, r//2, beta_bound, prec) beta_bound = refine_beta_bound(P, Q, E, mu, a + r//2, b + r//2, r//2, beta_bound, prec) - + return beta_bound end diff --git a/src/EllCrv/Isogeny.jl b/src/EllCrv/Isogeny.jl index c8e5ae7c41..f8e9d7e363 100644 --- a/src/EllCrv/Isogeny.jl +++ b/src/EllCrv/Isogeny.jl @@ -62,12 +62,12 @@ function is_kernel_polynomial(E::EllCrv{T}, f::PolyElem{T}, check::Bool = false) tor_2 = division_polynomial_univariate(E, 2)[1] ker_G = gcd(f, tor_2) - + #2-torsion can only be of degree 0, 1 or 3 if degree(ker_G) == 2 return false end - + #Kernel can't be cyclic if we have full 2-torsion if degree(ker_G)==3 && check return false @@ -82,13 +82,13 @@ function is_kernel_polynomial(E::EllCrv{T}, f::PolyElem{T}, check::Bool = false) tor_2 = division_polynomial_univariate(E, 2)[1] ker_G = gcd(f, tor_2) end - - #Now we check if the corresponding isogeny factors through some + + #Now we check if the corresponding isogeny factors through some #multiplication by m map. d = ZZRingElem(2*degree(f) + 1) n = Hecke.squarefree_part(d) m = sqrt(div(d, n)) - + M = divisors(m) c = length(M) psi_m = division_polynomial_univariate(E, M[c])[1] @@ -103,11 +103,11 @@ function is_kernel_polynomial(E::EllCrv{T}, f::PolyElem{T}, check::Bool = false) if degree(psi_m) != 0 && check return false end - + phi_m = Isogeny(E, psi_m) f = push_through_isogeny(phi_m, f) E = codomain(phi_m) - + #Now the remainder should give us a cyclic isogeny #Check if the remaining cyclic part actually defines an isogeny d = 2*degree(f)+1 @@ -119,18 +119,18 @@ function is_kernel_polynomial(E::EllCrv{T}, f::PolyElem{T}, check::Bool = false) if !is_prime_cyclic_kernel_polynomial(E, p, f) return false end - + phi = Isogeny(E, ker_G) f = push_through_isogeny(phi, f) E = codomain(phi) d = 2*degree(f)+1 L = sort(prime_divisors(d)) end - + return true - + end - + @@ -375,11 +375,11 @@ function dual_isogeny(psi::Isogeny) E = domain(psi) d = degree(psi) K = base_field(E) - + if K(d) == 0 error("Cannot compute duals of separable isogenies yet.") end - + psi_d = division_polynomial_univariate(E, d)[1] psinew = push_through_isogeny(psi, psi_d) psihat = isogeny_from_kernel(codomain(psi), psinew) @@ -387,16 +387,16 @@ function dual_isogeny(psi::Isogeny) trans = isomorphism(codomain(psihat), E) psihat_up_to_auto = psihat * trans - + #Compute the first coefficients of psi and psihat scalar_psi = coeff(formal_isogeny(psi, 5), 1) scalar_psihat = coeff(formal_isogeny(psihat_up_to_auto, 5), 1) - + #Their product needs to be equal to the degree of phi (i.e. the first coefficient of the multiplication by m map) scalar = scalar_psi*scalar_psihat//d if scalar != one(base_field(E)) aut_E = automorphism_list(E) - + for s in aut_E u = isomorphism_data(s)[4] if u == scalar @@ -417,15 +417,15 @@ Return the dual of frobenius. function dual_of_frobenius(E) supsing = is_supersingular(E) - + a1, a2, a3, a4, a6 = a_invars(E) f = Isogeny(E) K = base_field(E) p = characteristic(K) - + f.codomain = E f.degree = p - + p = characteristic(base_field(E)) phim = multiplication_by_m_map(E, p) psi = isogeny_map_psi(phim) @@ -434,34 +434,34 @@ function dual_of_frobenius(E) Kxy = parent(omega) y = gen(Kxy) x = gen(base_ring(Kxy)) - + #The omega part is written as omega0(x^p) + y^p*(omega1(x^p)) #We find omega1 by dividing the y-part by the y-part of y^p mod f(x, y) #Then we use this to compute omega0. - + #If the curve is supersingular the dual of Frobenius will be an automorphism composed with Frobenius, #so we write y^(p^2) = f(x) +y*g(x) to find the automorphism - #Otherwise it will be a separable isogeny and it suffices to find f and g such that + #Otherwise it will be a separable isogeny and it suffices to find f and g such that #y^p = = f(x) +y*g(x) if supsing yp = lift(residue_ring(Kxy, y^2 +a1*x*y + a3*y - x^3 - a2*x^2 - a4*x -a6)(y^(p^2))) else yp = lift(residue_ring(Kxy, y^2 +a1*x*y + a3*y - x^3 - a2*x^2 - a4*x -a6)(y^p)) end - + omega1 = divexact(coefficients(omega)[1], coefficients(yp)[1]) omega0 = coefficients(omega)[0] - coefficients(yp)[0] * omega1 - + dual_psi = defrobenify(psi, p) dual_phi = defrobenify(phi, p) - + dual_omega0= defrobenify(omega0, p) dual_omega1= defrobenify(omega1, p) - + f.coordx = dual_phi//(dual_psi)^2 f.coordy = (dual_omega0 + dual_omega1*y)//(Kxy(dual_psi))^3 - + f.psi = dual_psi f.header = MapHeader(E, f.codomain) if supsing @@ -526,7 +526,7 @@ function defrobenify(f::RingElem, p, rmax::Int = -1) if is_zero(f) return f end - + p = Int(p) nonzerocoeffs = [coefficients(f)[i]!=0 ? i : 0 for i in (0:p:degree(f))] pr = gcd(nonzerocoeffs) @@ -536,7 +536,7 @@ function defrobenify(f::RingElem, p, rmax::Int = -1) end r = valuation(pr, p) - if rmax >= 0 && rmax < r + if rmax >= 0 && rmax < r r = rmax pr = p^r end @@ -590,16 +590,16 @@ function rational_maps(I::Isogeny) fdenom = to_bivariate(denominator(I.coordx)) gnum = to_bivariate(numerator(I.coordy)) gdenom = to_bivariate(denominator(I.coordy)) - + Ix = fnum//fdenom Iy = gnum//gdenom - + return [Ix, Iy, one(parent(Ix))] end function Base.getindex(f::Isogeny, i::Int) @req 1 <= i <= 3 "Index must be 1, 2 or 3" - + return rational_maps(f)[i] end @@ -614,7 +614,7 @@ function show(io::IO, f::Isogeny) E1 = domain(f) E2 = codomain(f) fx, fy = rational_maps(f) - print(io, "Isogeny from + print(io, "Isogeny from $(E1) to \n $(E2) given by \n (x : y : 1) -> ($(fx) : $(fy) : 1 )") @@ -642,8 +642,8 @@ function compose(I1::Isogeny, I2::Isogeny) newx_overy = Rxy(numerator(newx))//Rxy(denominator(newx)) - - + + tempomega_num = numerator(I2.coordy) tempomega_denom = denominator(I2.coordy) @@ -696,7 +696,7 @@ function compose(fs::Vector{Isogeny}) end function ^(phi::Isogeny, n::Int) - + res = identity_isogeny(E) for i in (1:n) @@ -888,21 +888,21 @@ function to_bivariate(f::AbstractAlgebra.Generic.Poly{S}) where S<:PolyElem{T} w Kxy, (x, y) = polynomial_ring(R, ["x","y"]) cx = coefficients(f) - + newf = zero(Kxy) for i in (0:length(cx)) newf += cx[i](x)*y^i end - + return newf end function to_bivariate(f::PolyElem{T}) where T<:FieldElem K = base_ring(f) - + Kxy, (x, y) = polynomial_ring(K, ["x","y"]) - + return f(x) end diff --git a/src/EllCrv/Isomorphisms.jl b/src/EllCrv/Isomorphisms.jl index 88ea6e5a0c..5d8c4d73ed 100644 --- a/src/EllCrv/Isomorphisms.jl +++ b/src/EllCrv/Isomorphisms.jl @@ -599,7 +599,7 @@ end function Base.getindex(f::EllCrvIso, i::Int) @req 1 <= i <= 3 "Index must be 1, 2 or 3" - + return rational_maps(f)[i] end @@ -626,7 +626,7 @@ function show(io::IO, f::EllCrvIso) E1 = domain(f) E2 = codomain(f) fx, fy, fz = rational_maps(f) - print(io, "Isomorphism from + print(io, "Isomorphism from $(E1) to \n $(E2) given by \n (x : y : 1) -> ($(fx) : $(fy) : $(fz) )") diff --git a/src/EllCrv/LocalData.jl b/src/EllCrv/LocalData.jl index 521ff4d7bb..32ecfa27fa 100644 --- a/src/EllCrv/LocalData.jl +++ b/src/EllCrv/LocalData.jl @@ -811,19 +811,19 @@ struct KodairaSymbol if K[1]!='I' error("String does not represent a valid Kodaira symbol.") end - + n = lastindex(K) - + if K[n]=='*' m = parse(Int, K[2:n-1]) return KodairaSymbol(-4 - m) else m = parse(Int, K[2:n]) return KodairaSymbol(4 + m) - end - + end + error("String does not represent a valid Kodaira symbol.") - + end end @@ -879,17 +879,17 @@ function ==(K::KodairaSymbol, s::String) if s[1] != 'I' error("String does not represent a valid Kodaira symbol.") end - + n = lastindex(s) if s[n]=='*' m = parse(Int, s[2:n-1]) - return K.ksymbol == -4 - m + return K.ksymbol == -4 - m else m = parse(Int, s[2:n]) return K.ksymbol == 4 + m - end - + end + error("String does not represent a valid Kodaira symbol.") end @@ -897,7 +897,7 @@ end function show(io::IO, K::KodairaSymbol) m = K.ksymbol - if m == 1 + if m == 1 print(io, "I0") elseif m == -1 print(io, "I0*") @@ -914,15 +914,15 @@ function show(io::IO, K::KodairaSymbol) elseif m == -4 print(io, "IV*") end - - if m > 4 + + if m > 4 m = m - 4 print(io, "I$(m)") elseif m < -4 m = m + 4 m = -m print(io, "I$(m)*") - end + end end @@ -1026,7 +1026,7 @@ function reduction_type(E::EllCrv{QQFieldElem}, p) return "Good" end - if Kp.ksymbol > 4 + if Kp.ksymbol > 4 if split return "Split multiplicative" else @@ -1052,7 +1052,7 @@ function reduction_type(E::EllCrv{nf_elem}, p::NfOrdIdl) return "Good" end - if Kp.ksymbol > 4 + if Kp.ksymbol > 4 if split return "Split multiplicative" else diff --git a/src/EllCrv/MinimalModels.jl b/src/EllCrv/MinimalModels.jl index fa21390c5d..7734cbb9be 100644 --- a/src/EllCrv/MinimalModels.jl +++ b/src/EllCrv/MinimalModels.jl @@ -669,7 +669,7 @@ end # for k in range(n): # u*=U.gen(k).value()^min(e[k] for e in E) # return u -# +# # function _factor_nf(n::QQFieldElem) f = factor(ZZ, n) @@ -725,7 +725,7 @@ function _factor_rational_function_field(p) facpn = factor(pn) # pn is primitive and integral # we assume R is a UFD - + u = unit(facpn) for (p, e) in facpn uu, pp = _make_primitive(p) diff --git a/src/EllCrv/Misc.jl b/src/EllCrv/Misc.jl index 1d4bc92416..3e9ec59d00 100644 --- a/src/EllCrv/Misc.jl +++ b/src/EllCrv/Misc.jl @@ -119,9 +119,9 @@ function quadroots(a::nf_elem, b::nf_elem, c::nf_elem, pIdeal:: NfOrdIdl) R = order(pIdeal) F, phi = residue_field(R, pIdeal) P, x = polynomial_ring(F, "x", cached = false) - + t = [phi(R(numerator(s)))//phi(R(denominator(s))) for s in [a, b, c]] - + f = t[1]*x^2 + t[2]*x + t[3] if degree(f) == -1 @@ -179,7 +179,7 @@ function nrootscubic(b::nf_elem, c::nf_elem, d::nf_elem, pIdeal:: NfOrdIdl) R = order(pIdeal) F, phi = residue_field(R, pIdeal) P, x = polynomial_ring(F, "x", cached = false) - + t = [phi(R(numerator(s)))//phi(R(denominator(s))) for s in [b,c,d]] f = x^3 + t[1]*x^2 + t[2]*x + t[3] diff --git a/src/EllCrv/Models.jl b/src/EllCrv/Models.jl index 31f7647444..1da30dc30b 100644 --- a/src/EllCrv/Models.jl +++ b/src/EllCrv/Models.jl @@ -291,7 +291,7 @@ function is_integral_model(E::EllCrv{T}) where T<:Union{QQFieldElem, nf_elem} end @doc raw""" - is_local_integral_model(E::EllCrv{nf_elem}, P::NfOrdIdl) -> Bool + is_local_integral_model(E::EllCrv{nf_elem}, P::NfOrdIdl) -> Bool Given an elliptic curve $E$ over a number field $K$ and a prime ideal, return true if $E$ is a local integral model of $E$. diff --git a/src/EllCrv/Pairings.jl b/src/EllCrv/Pairings.jl index 5ae5b8ad84..931887e8d9 100644 --- a/src/EllCrv/Pairings.jl +++ b/src/EllCrv/Pairings.jl @@ -51,30 +51,30 @@ function straight_line(P::EllCrvPt{T}, Q::EllCrvPt{T}, R::EllCrvPt{T}) where T end E = parent(P) K = base_field(E) - + if is_infinite(P) || is_infinite(Q) if P == Q #Degenerate case return one(K) end - + if is_infinite(P) #Vertical line return R[1] - Q[1] end - + if is_infinite(Q) - #Vertical line + #Vertical line return R[1] - P[1] end end - + if P == Q #Line tangent to P a1, a2, a3, a4, a6 = a_invars(E) num = 3*P[1]^2 + 2*a2*P[1] + a4 - a1*P[2] denom = 2*P[2] + a1*P[1] + a3 - + if denom == 0 return R[1] - P[1] end @@ -88,7 +88,7 @@ function straight_line(P::EllCrvPt{T}, Q::EllCrvPt{T}, R::EllCrvPt{T}) where T slope = (Q[2] - P[2])// (Q[1] - P[1]) end end - + return R[2] - P[2] - slope * (R[1] - P[1]) end @@ -99,17 +99,17 @@ function _try_miller(P::EllCrvPt{T}, Q::EllCrvPt{T}, n::Int) where T @req parent(P) == parent(Q) "P and Q need to lie on the same curve" @req is_finite(Q) "Q must be finite" @req n != 0 "n must be non-zero" - + n_negative = n < 0 - if n_negative + if n_negative n = -n end - + t = one(base_field(parent(P))) V = P nbinary = digits(n, base=2) i = nbits(n) - 1 - + while i > 0 S = 2 * V ell = straight_line(V, V, Q) @@ -131,7 +131,7 @@ function _try_miller(P::EllCrvPt{T}, Q::EllCrvPt{T}, n::Int) where T end i = i - 1 end - + if n_negative vee = straight_line(V, -V, Q) if iszero(vee) || iszero(t) @@ -139,7 +139,7 @@ function _try_miller(P::EllCrvPt{T}, Q::EllCrvPt{T}, n::Int) where T end t = 1//(t*vee) end - return true, t + return true, t end @doc raw""" @@ -154,11 +154,11 @@ function weil_pairing(P::EllCrvPt{T}, Q::EllCrvPt{T}, n::Int) where T O = infinity(E) @req E == parent(Q) "P and Q need to be points on the same curve." @req n*P == O && n*Q == O "P and Q need to be n-torsion points." - + if P == Q || P == O || Q == O return one(K) end - + fl, m1 = _try_miller(P, Q, n) fl2, m2 = _try_miller(Q, P, n) if !fl || !fl2 @@ -173,7 +173,7 @@ end Given an $n$-torsion point $P$ and another point $Q$ on an elliptic curve over a finite field $K$. Return the Tate pairing $t_n(P, Q)$ by returning a -representative of the coset modulo $n$-th powers. +representative of the coset modulo $n$-th powers. See also [reduced_tate_pairing](@ref) if one wants an $n$-th root. """ @@ -211,12 +211,12 @@ See also [`tate_pairing`](@ref). function reduced_tate_pairing(P::EllCrvPt, Q::EllCrvPt, n) E = parent(P) @req E == parent(Q) "P and Q need to be points on the same curve." - K = base_field(E) + K = base_field(E) d = degree(K) q = order(K) @req divisible(q - 1, n) "K needs to contain the nth roots of unity." @req n*P == infinity(E) "P must be of order n." - + e = div(q - 1, n) return tate_pairing(P, Q, n)^e end diff --git a/src/EllCrv/RationalPointSearch.jl b/src/EllCrv/RationalPointSearch.jl index e607da4ddd..c3746a1722 100644 --- a/src/EllCrv/RationalPointSearch.jl +++ b/src/EllCrv/RationalPointSearch.jl @@ -11,8 +11,8 @@ #Algorithm based on ratpoints by Michael Stoll #Things Stoll does that we (may not) have yet (or may not want): -# - Sieving a second round. After the initial sieving round -# Stoll sieves the surviving chunks a second time using a different set of primes +# - Sieving a second round. After the initial sieving round +# Stoll sieves the surviving chunks a second time using a different set of primes # - Use gcd and Jacobi symbol routines that rely (almost) entirely on differences and # shifts. (Not sure if we have this already (or need this) # - Compute the residue classes of b modulo the various sieving primes by addition of the @@ -37,7 +37,7 @@ const _primes_for_sieve = find_points((coefficients::Vector{ZZRingElem}, bound::IntegerUnion) -> ArbField Given a list of coefficients [a_0, a_1, ..., a_n] and a bound, -return a list of points (x, y) satisfying +return a list of points (x, y) satisfying y^2 = a_0 + a_1*x + ... +a_n*x^n where writing x = a//b with gcd(a, b) = 1 the points are bounded by max(|a|, |b|) <= bound. @@ -82,7 +82,7 @@ function find_points(E::EllCrv{QQFieldElem}, bound::IntegerUnion, N = 2^14, P = end coeffs = map(numerator, coeffs) points = _find_points(coeffs, bound, N, P, Pfirst) - + if transform map!(t-> transform_ys(t, a1, a3), points, points) else @@ -108,7 +108,7 @@ function _find_points(coefficients::Vector, bound::Union{Integer, ZZRingElem}, N #N is size of chunks in which we subdivide #P is the number of primes we consider for sieving #Pfirst is the number of optimal primes we consider - + @req coefficients[end] != 0 "Leading coefficient needs to be non-zero" # Cache for some BitVector operations @@ -125,7 +125,7 @@ function _find_points(coefficients::Vector, bound::Union{Integer, ZZRingElem}, N primes = _primes_for_sieve[1:P] # Define the polynomial because we like to evaluate it - + # First decide whether to reverse the polynomial or not n = length(coefficients) @@ -143,20 +143,20 @@ function _find_points(coefficients::Vector, bound::Union{Integer, ZZRingElem}, N while tempcoeff[1] == 0 popfirst!(tempcoeff) end - + if isodd(length(tempcoeff)) && tempcoeff[1] < coefficients[n] reverse_polynomial = true coefficients = reverse!(tempcoeff) end end #If f is of even degree, reverse the polynomial if it would lead to better results - else + else #First compute the coefficients of the reciprocal polynomial potential_coefficients = reverse(coefficients) while potential_coefficients[end] == 0 pop!(potential_coefficients) end - + #If the reciprocal polynomial has odd degree we reverse if iseven(length(potential_coefficients)) reverse_polynomial = true @@ -165,29 +165,29 @@ function _find_points(coefficients::Vector, bound::Union{Integer, ZZRingElem}, N #Store the leading coefficients for later comparison old_lead = coefficients[end] new_lead = potential_coefficients[end] - end + end end - + g = Hecke.Globals.Qx(coefficients) - + n = length(coefficients) odd_degree = iseven(n) lead_coeff = coefficients[end] - + # Now first compute the primes we will use for sieving # We use those primes such that mod p there are few points - + # Take the Pfirst primes that are most optimal for sieving best_primes = Tuple{Int, QQFieldElem}[] exclude_denom= Int[] exclude_denom_new = Int[] - + for p in primes F = GF(p, cached = false) order = Hecke.order_via_exhaustive_search(map(F, coefficients)) push!(best_primes, (p, order//p)) - - if !odd_degree + + if !odd_degree if !is_square(F(old_lead)) push!(exclude_denom, p) end @@ -196,13 +196,13 @@ function _find_points(coefficients::Vector, bound::Union{Integer, ZZRingElem}, N end end end - - if !odd_degree - #We are going to exclude the denominators b divided by p for which the + + if !odd_degree + #We are going to exclude the denominators b divided by p for which the #squarefree part of the leading coefficient is a non-square mod p. #To maximize the number of denominators we use the euler_phi function as - #a heuristic. - + #a heuristic. + #Maybe simply compute the size of B for both potential leading coefficients? #This test doesn't include the effect of prime divisors of leading coefficient old_d = prod(exclude_denom;init = one(ZZRingElem)) @@ -213,7 +213,7 @@ function _find_points(coefficients::Vector, bound::Union{Integer, ZZRingElem}, N if old_score > new_score reverse_polynomial = true coefficients = potential_coefficients - + exclude_denom = exclude_denom_new end end @@ -222,9 +222,9 @@ function _find_points(coefficients::Vector, bound::Union{Integer, ZZRingElem}, N primes = Int[p for (p,q) in best_primes[1:Pfirst]] n = length(coefficients) - + lead_coeff = coefficients[n] - + #H[m][n] contains sieving info for the residue class k-1 mod m H = Vector{BitVector}[] H_temp = BitVector[] @@ -240,15 +240,15 @@ function _find_points(coefficients::Vector, bound::Union{Integer, ZZRingElem}, N push!(H_temp, copy(p_sieve[1])) push!(p_starts, mod(-bound, p)) end - + two_adic_info = mod16_check_arrays(coefficients) - + #candidates = fill(trues(N), H_parts) candidates = BitVector[trues(N) for i in 1:(H_parts + 1)] res = Tuple{QQFieldElem, QQFieldElem, QQFieldElem}[] - + #Determine the set of denumerators b #if the polynomial is odd we can restrict the possible denominators to the ones of the form a*b^2 @@ -263,7 +263,7 @@ function _find_points(coefficients::Vector, bound::Union{Integer, ZZRingElem}, N end B = collect(filter!(t -> t <= bound, BB)) else - #Exclude denominators divisible py p^k for which the leading coefficient is not a + #Exclude denominators divisible py p^k for which the leading coefficient is not a #square mod p^k. (k minimal in this way) for p in prime_divisors(lead_coeff) vp = valuation(lead_coeff, p) @@ -285,20 +285,20 @@ function _find_points(coefficients::Vector, bound::Union{Integer, ZZRingElem}, N end gcdb = gcd(b, 2*lead_coeff) - + while gcdb != 1 b = divexact(b, gcd(b, 2*lead_coeff)) gcdb = gcd(b, 2*lead_coeff) end - + if jacobi_symbol(lead_coeff, b) != 1 return false end - + return true end end - + neg_ints = negative_intervals(g) left = neg_ints[1] @@ -306,51 +306,51 @@ function _find_points(coefficients::Vector, bound::Union{Integer, ZZRingElem}, N right = neg_ints[3] interval_bounds = QQFieldElem[] - + if !isempty(left) push!(interval_bounds, max(-bound, left[1])) else push!(interval_bounds, -bound) end - + for I in intervals push!(interval_bounds, I[1]) push!(interval_bounds, I[2]) end - + if !isempty(right) push!(interval_bounds, min(bound, right[1])) else push!(interval_bounds, bound) end - + #Delete any intervals that are already out of bounds out_of_bounds = interval_bounds[2] < -bound while out_of_bounds - + popfirst!(interval_bounds) popfirst!(interval_bounds) out_of_bounds = interval_bounds[2] < -bound end - + out_of_bounds = interval_bounds[end - 1] > bound while out_of_bounds - + pop!(interval_bounds) pop!(interval_bounds) out_of_bounds = interval_bounds[end - 1] < bound end - + #Add point(s) at infinity of desingularized projective closure #If the degree of f is odd we have y^2 = a(n-1)*x^(n - 1)*z + ... + a0*z^n, so (1:0:0) is a point - if odd_degree_original + if odd_degree_original push!(res, (one(QQFieldElem), zero(QQFieldElem), zero(QQFieldElem))) #If the degree of f is even we have y^2 = an*x^n + ... + a0*z^n, so (1:1:0) and (1:-1:0) are points if and only if an is a square elseif is_square(leading_coefficient(f)) push!(res, (one(QQFieldElem), one(QQFieldElem), zero(QQFieldElem))) push!(res, (one(QQFieldElem), -one(QQFieldElem), zero(QQFieldElem))) end - + if reverse_polynomial points_with_x!(res, zero(QQFieldElem), f) end @@ -374,50 +374,50 @@ function _find_points_in_interval!(res, f, coefficients::Vector, primes, H_tripl else H = H_triple[case] end - + #Compute the start of the interval and the end of the interval, but make sure it is within height bounds start_interval = min(max(ceil(ZZRingElem, b*left_bound), -bound), bound) end_interval = max(min(floor(ZZRingElem, b*right_bound), bound), -bound) - - + + # If we only consider odd or even numerators if case > 1 - + #Make sure starting bit corresponds to even numerator if case = 2 and odd if case = 3 if isodd(start_interval + case) start_interval += 1 end - + #Range is divided by 2 when we only consider odd or even numerators numerator_range = ceil(Int, (1 - start_interval + end_interval)// 2) else numerator_range = 1 + - start_interval + end_interval end H_parts = Int(div(numerator_range, N)) - + rest = rem(numerator_range, N) - + for i in 1:H_parts fill!(candidates[i], true) end candidates[H_parts + 1] = BitVector(x <= rest for x in 1:N) - + for i in 1:length(primes) p = @inbounds primes[i] - + if p < N shift = -rem(N, p) else shift = -N end - + k = mod(b, p) if case == 1 offset = Int(mod(start_interval, p)) elseif case == 2 temp = Int(mod(start_interval, p)) if iseven(temp) - offset = divexact(temp, 2) + offset = divexact(temp, 2) else offset = div(p, 2) + 1 + div(temp, 2) end @@ -426,7 +426,7 @@ function _find_points_in_interval!(res, f, coefficients::Vector, primes, H_tripl if iseven(temp) offset = mod(div(p, 2) + divexact(temp, 2), p) else - offset = div(temp, 2) + offset = div(temp, 2) end end k = mod(b, p) @@ -434,7 +434,7 @@ function _find_points_in_interval!(res, f, coefficients::Vector, primes, H_tripl #Need to shift by 1 as H[i][k] corresponds to k-1 mod p p_sieve = @inbounds H[i][k + 1] p_sieve_temp = @inbounds H_temp[i] - + for j in 1:(H_parts + 1) c = candidates[j] @@ -460,14 +460,14 @@ function _find_points_in_interval!(res, f, coefficients::Vector, primes, H_tripl for j in 1:cnt a = _a + I[j] if gcd(a, b) == 1 - if reverse_polynomial - if a != 0 + if reverse_polynomial + if a != 0 x = QQFieldElem(b//a) else continue end else - + x = QQFieldElem(a//b) end points_with_x!(res, x, f) @@ -487,8 +487,8 @@ function _find_points_in_interval!(res, f, coefficients::Vector, primes, H_tripl for j in 1:cnt a = _a + 2*(I[j] - 1) if gcd(a, b) == 1 - if reverse_polynomial - if a != 0 + if reverse_polynomial + if a != 0 x = QQFieldElem(b//a) else continue @@ -518,24 +518,24 @@ function prime_check_arrays(coeff::Vector{<: IntegerUnion}, p::Int, N) p_part = Vector{Vector{Bool}}(undef, p) p_part_odd = Vector{Vector{Bool}}(undef, p) p_part_even = Vector{Vector{Bool}}(undef, p) - - + + az = Vector{elem_type(F)}(undef, n + 1) _chunk = Vector{Bool}(undef, length(F)) for t in (0:p - 1) z = F(t) zpowers = powers(z, n) - #a[i+1] corresponds to a_i above + #a[i+1] corresponds to a_i above for i in 0:n az[i + 1] = a[i + 1] * zpowers[n - i + 1] end - + chunk = _chunk #for (j, x) in enumerate(F) # @inbounds chunk[j] = issquare(sum([az[i + 1]*x^i for i in (0:n)])) #end chunk = Bool[issquare(sum([az[i + 1]*x^i for i in (0:n)])) for x in F] - chunk_odd = vcat(chunk[2:2:p], chunk[1:2:p]) + chunk_odd = vcat(chunk[2:2:p], chunk[1:2:p]) chunk_even = vcat(chunk[1:2:p], chunk[2:2:p]) #Pad the BitArray to have chunks that are at least big enough to do a broadcasted & with @@ -582,7 +582,7 @@ function prime_check_arrays_2(coeff::Vector{<: IntegerUnion}, p::Int, N, C) p_part = Vector{BitVector}(undef, p) p_part_odd = Vector{BitVector}(undef, p) p_part_even = Vector{BitVector}(undef, p) - + temp = BitVector(undef, p) if p < N @@ -594,7 +594,7 @@ function prime_check_arrays_2(coeff::Vector{<: IntegerUnion}, p::Int, N, C) az = Vector{elem_type(F)}(undef, n + 1) collF = collect(F) - + temp = BitVector(undef, p) temp2 = BitVector(undef, p) temp3 = BitVector(undef, p) @@ -616,12 +616,12 @@ function prime_check_arrays_2(coeff::Vector{<: IntegerUnion}, p::Int, N, C) # compute z,...,z^{n+1} if z is even and # 1,...,z^n if n is even zpowers = _powers(z, isoddn, n) - + #a[i+1] corresponds to a_i above for i in 0:n az[i + 1] = a[i + 1] * zpowers[n - i + 1] end - + for i in 1:p temp4[i] = squares[_sum_up(az, collF[i], n).data + 1] end @@ -698,19 +698,19 @@ function mod16_check_arrays(coefficients::Vector{<: IntegerUnion}) a = map(R, coefficients) part_16 = Array{Int}(undef, 16) - + if isodd(n) d = 1 else d = 0 end - + # t odd for t in (1:2:15) z = R(t) - + #Projective closure - + #a[i+1] corresponds to a_i above chunk = BitArray(sum([a[i + 1]*x^i*z^(n - i + d) for i in (0:n)]) in map(R, [0,1,4,9]) for x in R) if chunk == falses(16) @@ -718,7 +718,7 @@ function mod16_check_arrays(coefficients::Vector{<: IntegerUnion}) else evens = [chunk[i] for i in (1:2:16)] odds = [chunk[i] for i in (2:2:16)] - + #Only even solutions if odds == falses(8) part_16[t+1] = 2 @@ -731,7 +731,7 @@ function mod16_check_arrays(coefficients::Vector{<: IntegerUnion}) end end end - + for t in (0:2:15) z = R(t) #a[i+1] corresponds to a_i above diff --git a/src/EllCrv/Torsion.jl b/src/EllCrv/Torsion.jl index d124f22623..6ff0179dd0 100644 --- a/src/EllCrv/Torsion.jl +++ b/src/EllCrv/Torsion.jl @@ -349,7 +349,7 @@ end torsion_bound(E::EllCrv{nf_elem}, n::Int) -> ZZRingElem -Bound the order of the torsion subgroup of $E by considering +Bound the order of the torsion subgroup of $E by considering the order of the reduction of $E$ modulo $n$ distinct primes with good reduction """ @@ -373,7 +373,7 @@ function torsion_bound(E::EllCrv{nf_elem}, n::Int) end end p = next_prime(p) - end + end return(ZZRingElem(bound)) end diff --git a/test/setup.jl b/test/setup.jl index 1f332f6c6d..981f4c856a 100644 --- a/test/setup.jl +++ b/test/setup.jl @@ -34,7 +34,7 @@ else end end -if (isdefined(Hecke, :_with_gap) && Hecke._with_gap) || +if (isdefined(Hecke, :_with_gap) && Hecke._with_gap) || (isdefined(Main, :_with_gap) && Main._with_gap) macro with_gap(ex) ex @@ -45,7 +45,7 @@ else end end -if (isdefined(Hecke, :_with_polymake) && Hecke._with_polymake) || +if (isdefined(Hecke, :_with_polymake) && Hecke._with_polymake) || (isdefined(Main, :_with_polymake) && Main._with_polymake) macro with_polymake(ex) ex