diff --git a/mpmath/functions/elliptic.py b/mpmath/functions/elliptic.py index 04d5534c..a78d2df3 100644 --- a/mpmath/functions/elliptic.py +++ b/mpmath/functions/elliptic.py @@ -1718,30 +1718,34 @@ def g2g3from(ctx, q=None, m=None, k=None, tau=None, qbar=None, "g2g3from", q, m, k, tau, qbar, g2, g3, omega1, omega2) with ctx.extraprec(10): if g2 is not None: - return +ctx.convert(g2), +ctx.convert(g3) - if omega1 is None: - tau = ctx.taufrom(q=q, m=m, k=k, tau=tau, qbar=qbar) - if ctx.im(tau) <= 0: - raise ValueError("g2g3from: tau must be in upper half-plane") - omega1 = ctx.one/2 - omega2 = tau/2 + g2 = ctx.convert(g2) + g3 = ctx.convert(g3) else: - omega1 = ctx.convert(omega1) - omega2 = ctx.convert(omega2) - if ctx.im(omega2/omega1) <= 0: - raise ValueError("g2g3from: omega ratio must be " - "in upper half-plane") - tau = omega2 / omega1 - q = ctx.qfrom(tau=tau) - j2 = ctx.jtheta(2, 0, q) - j3 = ctx.jtheta(3, 0, q) - factor = ctx.pi / (2 * omega1) - g2 = (ctx.mpf(4)/3) * factor**4 * (j2**8 - (j2*j3)**4 + j3**8) - g3 = ((ctx.mpf(8)/27) * factor**6 * - (j2**12 - (ctx.mpf(3)/2*j2**8*j3**4 + - ctx.mpf(3)/2*j2**4*j3**8) + - j3**12)) - return +g2, +g3 + if omega1 is None: + tau = ctx.taufrom(q=q, m=m, k=k, tau=tau, qbar=qbar) + if ctx.im(tau) <= 0: + raise ValueError("g2g3from: tau must be in upper " + "half-plane") + omega1 = ctx.one/2 + omega2 = tau/2 + else: + omega1 = ctx.convert(omega1) + omega2 = ctx.convert(omega2) + if ctx.im(omega2/omega1) <= 0: + raise ValueError("g2g3from: omega ratio must be " + "in upper half-plane") + tau = omega2 / omega1 + q = ctx.qfrom(tau=tau) + j2 = ctx.jtheta(2, 0, q) + j3 = ctx.jtheta(3, 0, q) + factor = ctx.pi / (2 * omega1) + g2 = ((ctx.mpf(4)/3) * factor**4 * + (j2**8 - (j2*j3)**4 + j3**8)) + g3 = ((ctx.mpf(8)/27) * factor**6 * + (j2**12 - (ctx.mpf(3)/2*j2**8*j3**4 + + ctx.mpf(3)/2*j2**4*j3**8) + + j3**12)) + return +g2, +g3 @defun def omega1omega2from(ctx, q=None, m=None, k=None, tau=None, qbar=None, @@ -1778,23 +1782,28 @@ def omega1omega2from(ctx, q=None, m=None, k=None, tau=None, qbar=None, """ _validate_weierstrass_parameter_args( "omega1omega2from", q, m, k, tau, qbar, g2, g3, omega1, omega2) - with ctx.extraprec(10): - if omega1 is not None: + if omega1 is not None: + with ctx.extraprec(10): omega1 = ctx.convert(omega1) omega2 = ctx.convert(omega2) if ctx.im(omega2/omega1) <= 0: raise ValueError("omega1omega2from: omega ratio must be " "in upper half-plane") - return +omega1, +omega2 - if g2 is None: + return +omega1, +omega2 + if g2 is None: + with ctx.extraprec(10): tau = ctx.taufrom(q=q, m=m, k=k, tau=tau, qbar=qbar) if ctx.im(tau) <= 0: raise ValueError("omega1omega2from: tau must be in upper " "half-plane") - return +(ctx.one/2), +(tau/2) + omega1 = ctx.one/2 + omega2 = tau/2 + return +omega1, +omega2 + with ctx.extraprec(10): g2 = ctx.convert(g2) g3 = ctx.convert(g3) + periods = None if g2 == 0: omegaA = (g3 ** (ctx.mpf(-1)/ctx.mpf(6)) * @@ -1850,41 +1859,52 @@ def omega1omega2from(ctx, q=None, m=None, k=None, tau=None, qbar=None, # For m <= 1/2 the standard rectangular basis is already # reduced. For m > 1/2 apply its S-transform directly. if m <= ctx.one/2: - return +real_period, +(ctx.j*imaginary_period) - return +(-ctx.j*imaginary_period), +real_period + periods = (real_period, + ctx.j*imaginary_period) + else: + periods = (-ctx.j*imaginary_period, real_period) + else: + # One real root and a conjugate pair. Real Cardano radicals + # followed by a quadratic transformation express both + # periods using complete elliptic integrals with real + # parameters. + root_discriminant = ctx.sqrt(-real_discriminant/1728) + exponent = ctx.mpf(1)/3 + u3 = real_g3/8 + root_discriminant + v3 = real_g3/8 - root_discriminant + u = ctx.sign(u3)*abs(u3)**exponent + v = ctx.sign(v3)*abs(v3)**exponent + real_root = u+v + root_real_part = 3*real_root/2 + root_imaginary_part = ctx.sqrt(3)*(u-v)/2 + H = ctx.sqrt(root_real_part**2 + + root_imaginary_part**2) + m = (H-root_real_part)/(2*H) + sqrt_H = ctx.sqrt(H) + real_period = ctx.ellipk(m)/sqrt_H + imaginary_part = ctx.ellipk(1-m)/(2*sqrt_H) - # One real root and a conjugate pair. Real Cardano radicals - # followed by a quadratic transformation express both periods - # using complete elliptic integrals with real parameters. - root_discriminant = ctx.sqrt(-real_discriminant/1728) - exponent = ctx.mpf(1)/3 - u3 = real_g3/8 + root_discriminant - v3 = real_g3/8 - root_discriminant - u = ctx.sign(u3)*abs(u3)**exponent - v = ctx.sign(v3)*abs(v3)**exponent - real_root = u+v - root_real_part = 3*real_root/2 - root_imaginary_part = ctx.sqrt(3)*(u-v)/2 - H = ctx.sqrt(root_real_part**2 + root_imaginary_part**2) - m = (H-root_real_part)/(2*H) - sqrt_H = ctx.sqrt(H) - real_period = ctx.ellipk(m)/sqrt_H - imaginary_part = ctx.ellipk(1-m)/(2*sqrt_H) - - # These bases directly implement the documented fundamental- - # domain and simultaneous-sign conventions in each real - # symmetry region, so no generic PSL(2,Z) reduction is needed. - if real_g2 > 0: - if real_g3 > 0: - return (+real_period, - +(real_period/2 + ctx.j*imaginary_part)) - return (+(-2*ctx.j*imaginary_part), - +(real_period/2-ctx.j*imaginary_part)) - if real_g3 > 0: - return (+(real_period/2+ctx.j*imaginary_part), - +(-real_period/2+ctx.j*imaginary_part)) - return (+(real_period/2-ctx.j*imaginary_part), - +(real_period/2+ctx.j*imaginary_part)) + # These bases directly implement the documented + # fundamental-domain and simultaneous-sign conventions in + # each real symmetry region, so no generic PSL(2,Z) + # reduction is needed. + if real_g2 > 0: + if real_g3 > 0: + periods = ( + real_period, + real_period/2 + ctx.j*imaginary_part) + else: + periods = ( + -2*ctx.j*imaginary_part, + real_period/2-ctx.j*imaginary_part) + elif real_g3 > 0: + periods = ( + real_period/2+ctx.j*imaginary_part, + -real_period/2+ctx.j*imaginary_part) + else: + periods = ( + real_period/2-ctx.j*imaginary_part, + real_period/2+ctx.j*imaginary_part) else: # Solve the original cubic directly. This obtains both the # elliptic parameter and its scale together, avoiding inverse @@ -1907,35 +1927,39 @@ def omega1omega2from(ctx, q=None, m=None, k=None, tau=None, qbar=None, (2*ctx.agm(1, ctx.sqrt(m))*sqrt_D)) tau = omegaB/omegaA - if g2 != 0 and g3 != 0: - a, b, c, d = ctx._reduce_psl2z(tau) - omega1 = d*omegaA + c*omegaB - omega2 = b*omegaA + a*omegaB + if periods is None: + if g2 != 0 and g3 != 0: + a, b, c, d = ctx._reduce_psl2z(tau) + omega1 = d*omegaA + c*omegaB + omega2 = b*omegaA + a*omegaB + else: + omega1, omega2 = omegaA, omegaB + + # Remove sub-precision components introduced by inverse-j branch + # arithmetic before applying conventions on symmetry boundaries. + omega1 = ctx.chop(omega1) + omega2 = ctx.chop(omega2) + tau = omega2/omega1 + # Identify equivalent points on the vertical and circular + # boundaries of the fundamental domain without recomputing the + # period ratio. + if ctx.almosteq(ctx.re(tau), -ctx.one/2): + omega2 += omega1 + tau += 1 + if ctx.almosteq(abs(tau), ctx.one) and ctx.re(tau) < 0: + omega1, omega2 = omega2, -omega1 + + # The invariants do not distinguish a basis from its simultaneous + # negation. Place omega1 in the right half-plane, with the negative + # imaginary axis included as its boundary. + real = ctx.re(omega1) + if ((real < 0 and not ctx.almosteq(real, 0)) or + (ctx.almosteq(real, 0) and ctx.im(omega1) > 0)): + omega1 = -omega1 + omega2 = -omega2 else: - omega1, omega2 = omegaA, omegaB - - # Remove sub-precision components introduced by inverse-j branch - # arithmetic before applying conventions on symmetry boundaries. - omega1 = ctx.chop(omega1) - omega2 = ctx.chop(omega2) - tau = omega2/omega1 - # Identify equivalent points on the vertical and circular boundaries - # of the fundamental domain without recomputing the period ratio. - if ctx.almosteq(ctx.re(tau), -ctx.one/2): - omega2 += omega1 - tau += 1 - if ctx.almosteq(abs(tau), ctx.one) and ctx.re(tau) < 0: - omega1, omega2 = omega2, -omega1 - - # The invariants do not distinguish a basis from its simultaneous - # negation. Place omega1 in the right half-plane, with the negative - # imaginary axis included as its boundary. - real = ctx.re(omega1) - if ((real < 0 and not ctx.almosteq(real, 0)) or - (ctx.almosteq(real, 0) and ctx.im(omega1) > 0)): - omega1 = -omega1 - omega2 = -omega2 - return +omega1, +omega2 + omega1, omega2 = periods + return +omega1, +omega2 # ============================================================================ diff --git a/mpmath/tests/test_elliptic.py b/mpmath/tests/test_elliptic.py index fd2af26d..b013241a 100644 --- a/mpmath/tests/test_elliptic.py +++ b/mpmath/tests/test_elliptic.py @@ -1065,6 +1065,32 @@ def test_weierstrass_half_periods_direct_agm_roundtrip(): assert abs(tau) >= 1 - tol +def test_weierstrass_parameter_conversion_normalization(): + with mp.workprec(100): + high_precision = sqrt(2) + + with mp.workprec(53): + invariant_results = [ + g2g3from(g2=high_precision, g3=high_precision + 1), + g2g3from(tau=mpc('0.3', '1.2')), + g2g3from(omega1=high_precision, + omega2=j*high_precision), + ] + period_results = [ + omega1omega2from(omega1=high_precision, + omega2=j*high_precision), + omega1omega2from(tau=mpc('0.3', '1.2')), + omega1omega2from(g2=12, g3=1), + omega1omega2from(g2=60, g3=140), + omega1omega2from(g2=1 + 2*j, g3=3 - 4*j), + omega1omega2from(g2=1, g3=0), + omega1omega2from(g2=0, g3=1), + ] + + for result in invariant_results + period_results: + assert result == tuple(+value for value in result) + + def test_weierstrass_period_method_switch_is_continuous(): mp.dps = 50 threshold = ldexp(1, -20)