Merge pull request #1117 from HeskethGD/enhance-kleinj-tau
Add kleinjinv and support g and omega in taufrom
This commit is contained in:
@@ -63,3 +63,4 @@ Modular functions
|
||||
|
||||
.. autofunction:: mpmath.eta
|
||||
.. autofunction:: mpmath.kleinj
|
||||
.. autofunction:: mpmath.kleinjinv
|
||||
|
||||
@@ -46,6 +46,7 @@ qbarfrom = mp.qbarfrom
|
||||
ellipfun = mp.ellipfun
|
||||
jtheta = mp.jtheta
|
||||
kleinj = mp.kleinj
|
||||
kleinjinv = mp.kleinjinv
|
||||
eta = mp.eta
|
||||
|
||||
# Weierstrass elliptic functions
|
||||
|
||||
+176
-64
@@ -10,10 +10,11 @@ a function which satisfies
|
||||
f(z + 2 \omega_1) = f(z + 2 \omega_2) = f(z)
|
||||
|
||||
for some half-periods `\omega_1, \omega_2` with
|
||||
`\mathrm{Im}[\omega_1 / \omega_2] > 0`. The canonical elliptic
|
||||
functions are the Jacobi elliptic functions. More broadly, this section
|
||||
includes quasi-doubly periodic functions (such as the Jacobi theta
|
||||
functions) and other functions useful in the study of elliptic functions.
|
||||
`\mathrm{Im}[\omega_2 / \omega_1] > 0`. The classical examples include
|
||||
the Jacobi elliptic functions and the Weierstrass elliptic functions.
|
||||
More broadly, this section includes quasi-doubly periodic functions
|
||||
(such as the Jacobi theta functions) and other functions useful in the
|
||||
study of elliptic functions.
|
||||
|
||||
Many different conventions for the arguments of
|
||||
elliptic functions are in use. It is even standard to use
|
||||
@@ -24,6 +25,9 @@ must satisfy `|q| < 1`; the elliptic parameter `m` (an arbitrary
|
||||
complex number); the elliptic modulus `k` (an arbitrary complex
|
||||
number); and the half-period ratio `\tau`, which usually must
|
||||
satisfy `\mathrm{Im}[\tau] > 0`.
|
||||
Elliptic functions may also be parameterized by half-periods
|
||||
`\omega_1, \omega_2`, where `\tau = \omega_2 / \omega_1`, or by
|
||||
Weierstrass invariants `g_2, g_3`.
|
||||
These quantities can be expressed in terms of each other
|
||||
using the following relations:
|
||||
|
||||
@@ -50,8 +54,9 @@ number theory, which we here denote by q-bar:
|
||||
|
||||
\bar{q} = q^2 = e^{2 i \pi \tau}
|
||||
|
||||
For convenience, mpmath provides functions to convert
|
||||
between the various parameters (:func:`~mpmath.qfrom`, :func:`~mpmath.mfrom`,
|
||||
For convenience, mpmath provides functions to convert between the various
|
||||
parameters, including half-periods and Weierstrass invariants where
|
||||
applicable (:func:`~mpmath.qfrom`, :func:`~mpmath.mfrom`,
|
||||
:func:`~mpmath.kfrom`, :func:`~mpmath.taufrom`, :func:`~mpmath.qbarfrom`).
|
||||
|
||||
**References**
|
||||
@@ -120,10 +125,19 @@ def nome(ctx, m):
|
||||
v = ctx.mpc(0, v.imag)
|
||||
return v
|
||||
|
||||
def _validate_weierstrass_pair_args(funcname, g2, g3, omega1, omega2):
|
||||
if (g2 is None) != (g3 is None):
|
||||
raise ValueError("%s: must provide both g2 and g3" % funcname)
|
||||
if (omega1 is None) != (omega2 is None):
|
||||
raise ValueError("%s: must provide both omega1 and omega2" % funcname)
|
||||
|
||||
@defun_wrapped
|
||||
def qfrom(ctx, q=None, m=None, k=None, tau=None, qbar=None):
|
||||
def qfrom(ctx, q=None, m=None, k=None, tau=None, qbar=None,
|
||||
g2=None, g3=None, omega1=None, omega2=None):
|
||||
r"""
|
||||
Returns the elliptic nome `q`, given any of `q, m, k, \tau, \bar{q}`::
|
||||
Returns the elliptic nome `q`, given any of `q, m, k, \tau, \bar{q}`.
|
||||
It may also be computed from both Weierstrass invariants `g_2, g_3`,
|
||||
or both half-periods `\omega_1, \omega_2`::
|
||||
|
||||
>>> from mpmath import mp, qfrom, mfrom, kfrom, taufrom, qbarfrom
|
||||
>>> mp.dps = 25
|
||||
@@ -138,8 +152,11 @@ def qfrom(ctx, q=None, m=None, k=None, tau=None, qbar=None):
|
||||
(0.25 + 0.0j)
|
||||
>>> qfrom(qbar=qbarfrom(q=0.25))
|
||||
0.25
|
||||
>>> qfrom(omega1=1, omega2=0.5j)
|
||||
(0.2078795763507619085469556 + 0.0j)
|
||||
|
||||
"""
|
||||
_validate_weierstrass_pair_args("qfrom", g2, g3, omega1, omega2)
|
||||
if q is not None:
|
||||
return ctx.convert(q)
|
||||
if m is not None:
|
||||
@@ -150,12 +167,18 @@ def qfrom(ctx, q=None, m=None, k=None, tau=None, qbar=None):
|
||||
return ctx.expjpi(tau)
|
||||
if qbar is not None:
|
||||
return ctx.sqrt(qbar)
|
||||
if g2 is not None:
|
||||
return ctx.qfrom(tau=ctx.taufrom(g2=g2, g3=g3))
|
||||
if omega1 is not None:
|
||||
return ctx.qfrom(tau=ctx.taufrom(omega1=omega1, omega2=omega2))
|
||||
|
||||
@defun_wrapped
|
||||
def qbarfrom(ctx, q=None, m=None, k=None, tau=None, qbar=None):
|
||||
def qbarfrom(ctx, q=None, m=None, k=None, tau=None, qbar=None,
|
||||
g2=None, g3=None, omega1=None, omega2=None):
|
||||
r"""
|
||||
Returns the number-theoretic nome `\bar q`, given any of
|
||||
`q, m, k, \tau, \bar{q}`::
|
||||
`q, m, k, \tau, \bar{q}`, both Weierstrass invariants `g_2, g_3`,
|
||||
or both half-periods `\omega_1, \omega_2`::
|
||||
|
||||
>>> from mpmath import (mp, qbarfrom, qfrom, extraprec, mfrom,
|
||||
... kfrom, taufrom)
|
||||
@@ -171,8 +194,11 @@ def qbarfrom(ctx, q=None, m=None, k=None, tau=None, qbar=None):
|
||||
0.25
|
||||
>>> qbarfrom(tau=taufrom(qbar=0.25))
|
||||
(0.25 + 0.0j)
|
||||
>>> qbarfrom(omega1=1, omega2=0.5j)
|
||||
(0.04321391826377224977441774 + 0.0j)
|
||||
|
||||
"""
|
||||
_validate_weierstrass_pair_args("qbarfrom", g2, g3, omega1, omega2)
|
||||
if qbar is not None:
|
||||
return ctx.convert(qbar)
|
||||
if q is not None:
|
||||
@@ -183,12 +209,18 @@ def qbarfrom(ctx, q=None, m=None, k=None, tau=None, qbar=None):
|
||||
return nome(ctx, ctx.convert(k)**2) ** 2
|
||||
if tau is not None:
|
||||
return ctx.expjpi(2*tau)
|
||||
if g2 is not None:
|
||||
return ctx.qbarfrom(tau=ctx.taufrom(g2=g2, g3=g3))
|
||||
if omega1 is not None:
|
||||
return ctx.qbarfrom(tau=ctx.taufrom(omega1=omega1, omega2=omega2))
|
||||
|
||||
@defun_wrapped
|
||||
def taufrom(ctx, q=None, m=None, k=None, tau=None, qbar=None):
|
||||
def taufrom(ctx, q=None, m=None, k=None, tau=None, qbar=None,
|
||||
g2=None, g3=None, omega1=None, omega2=None):
|
||||
r"""
|
||||
Returns the elliptic half-period ratio `\tau`, given any of
|
||||
`q, m, k, \tau, \bar{q}`::
|
||||
`q, m, k, \tau, \bar{q}`, both Weierstrass invariants
|
||||
`g_2, g_3`, or both half-periods `\omega_1, \omega_2`::
|
||||
|
||||
>>> from mpmath import mp, taufrom, qfrom, mfrom, kfrom, qbarfrom
|
||||
>>> mp.dps = 25
|
||||
@@ -203,8 +235,11 @@ def taufrom(ctx, q=None, m=None, k=None, tau=None, qbar=None):
|
||||
(0.0 + 0.5j)
|
||||
>>> taufrom(qbar=qbarfrom(tau=0.5j))
|
||||
(0.0 + 0.5j)
|
||||
>>> taufrom(omega1=1, omega2=0.5j)
|
||||
(0.0 + 0.5j)
|
||||
|
||||
"""
|
||||
_validate_weierstrass_pair_args("taufrom", g2, g3, omega1, omega2)
|
||||
if tau is not None:
|
||||
return ctx.convert(tau)
|
||||
if m is not None:
|
||||
@@ -218,12 +253,25 @@ def taufrom(ctx, q=None, m=None, k=None, tau=None, qbar=None):
|
||||
if qbar is not None:
|
||||
qbar = ctx.convert(qbar)
|
||||
return ctx.log(qbar) / (2*ctx.pi*ctx.j)
|
||||
if g2 is not None:
|
||||
J = ctx.kleinj(g2=g2, g3=g3)
|
||||
return ctx.kleinjinv(J)
|
||||
if omega1 is not None:
|
||||
omega1 = ctx.convert(omega1)
|
||||
omega2 = ctx.convert(omega2)
|
||||
tau = omega2 / omega1
|
||||
if ctx.im(tau) <= 0:
|
||||
raise ValueError("taufrom: omega ratio must be in upper "
|
||||
"half-plane")
|
||||
return tau
|
||||
|
||||
@defun_wrapped
|
||||
def kfrom(ctx, q=None, m=None, k=None, tau=None, qbar=None):
|
||||
def kfrom(ctx, q=None, m=None, k=None, tau=None, qbar=None,
|
||||
g2=None, g3=None, omega1=None, omega2=None):
|
||||
r"""
|
||||
Returns the elliptic modulus `k`, given any of
|
||||
`q, m, k, \tau, \bar{q}`::
|
||||
`q, m, k, \tau, \bar{q}`, both Weierstrass invariants `g_2, g_3`,
|
||||
or both half-periods `\omega_1, \omega_2`::
|
||||
|
||||
>>> from mpmath import mp, kfrom, mfrom, qfrom, taufrom, qbarfrom
|
||||
>>> mp.dps = 25
|
||||
@@ -238,6 +286,8 @@ def kfrom(ctx, q=None, m=None, k=None, tau=None, qbar=None):
|
||||
(0.25 + 0.0j)
|
||||
>>> kfrom(qbar=qbarfrom(k=0.25))
|
||||
0.25
|
||||
>>> kfrom(omega1=1, omega2=0.5j)
|
||||
(0.985171431009416038689502 + 0.0j)
|
||||
|
||||
As `q \to 1` and `q \to -1`, `k` rapidly approaches
|
||||
`1` and `i \infty` respectively::
|
||||
@@ -251,6 +301,7 @@ def kfrom(ctx, q=None, m=None, k=None, tau=None, qbar=None):
|
||||
>>> kfrom(q=-1)
|
||||
(0.0 + infj)
|
||||
"""
|
||||
_validate_weierstrass_pair_args("kfrom", g2, g3, omega1, omega2)
|
||||
if k is not None:
|
||||
return ctx.convert(k)
|
||||
if m is not None:
|
||||
@@ -259,6 +310,11 @@ def kfrom(ctx, q=None, m=None, k=None, tau=None, qbar=None):
|
||||
q = ctx.expjpi(tau)
|
||||
if qbar is not None:
|
||||
q = ctx.sqrt(qbar)
|
||||
if q is None:
|
||||
if g2 is not None:
|
||||
return ctx.kfrom(tau=ctx.taufrom(g2=g2, g3=g3))
|
||||
if omega1 is not None:
|
||||
return ctx.kfrom(tau=ctx.taufrom(omega1=omega1, omega2=omega2))
|
||||
if q == 1:
|
||||
return q
|
||||
if q == -1:
|
||||
@@ -266,10 +322,12 @@ def kfrom(ctx, q=None, m=None, k=None, tau=None, qbar=None):
|
||||
return (ctx.jtheta(2,0,q)/ctx.jtheta(3,0,q))**2
|
||||
|
||||
@defun_wrapped
|
||||
def mfrom(ctx, q=None, m=None, k=None, tau=None, qbar=None):
|
||||
def mfrom(ctx, q=None, m=None, k=None, tau=None, qbar=None,
|
||||
g2=None, g3=None, omega1=None, omega2=None):
|
||||
r"""
|
||||
Returns the elliptic parameter `m`, given any of
|
||||
`q, m, k, \tau, \bar{q}`::
|
||||
`q, m, k, \tau, \bar{q}`, both Weierstrass invariants `g_2, g_3`,
|
||||
or both half-periods `\omega_1, \omega_2`::
|
||||
|
||||
>>> from mpmath import mp, mfrom, qfrom, kfrom, taufrom, qbarfrom, taylor
|
||||
>>> mp.dps = 25
|
||||
@@ -284,6 +342,8 @@ def mfrom(ctx, q=None, m=None, k=None, tau=None, qbar=None):
|
||||
(0.25 + 0.0j)
|
||||
>>> mfrom(qbar=qbarfrom(m=0.25))
|
||||
0.25
|
||||
>>> mfrom(omega1=1, omega2=0.5j)
|
||||
(0.9705627484771405856202647 + 0.0j)
|
||||
|
||||
As `q \to 1` and `q \to -1`, `m` rapidly approaches
|
||||
`1` and `-\infty` respectively::
|
||||
@@ -304,6 +364,7 @@ def mfrom(ctx, q=None, m=None, k=None, tau=None, qbar=None):
|
||||
[0.0, 16.0, -128.0, 704.0, -3072.0, 11488.0, -38400.0, 117632.0]
|
||||
|
||||
"""
|
||||
_validate_weierstrass_pair_args("mfrom", g2, g3, omega1, omega2)
|
||||
if m is not None:
|
||||
return m
|
||||
if k is not None:
|
||||
@@ -312,6 +373,11 @@ def mfrom(ctx, q=None, m=None, k=None, tau=None, qbar=None):
|
||||
q = ctx.expjpi(tau)
|
||||
if qbar is not None:
|
||||
q = ctx.sqrt(qbar)
|
||||
if q is None:
|
||||
if g2 is not None:
|
||||
return ctx.mfrom(tau=ctx.taufrom(g2=g2, g3=g3))
|
||||
if omega1 is not None:
|
||||
return ctx.mfrom(tau=ctx.taufrom(omega1=omega1, omega2=omega2))
|
||||
if q == 1:
|
||||
return ctx.convert(q)
|
||||
if q == -1:
|
||||
@@ -379,7 +445,8 @@ def ellipfun(ctx, kind, u=None, m=None, q=None, k=None, tau=None):
|
||||
return +v
|
||||
|
||||
@defun_wrapped
|
||||
def kleinj(ctx, tau=None, **kwargs):
|
||||
def kleinj(ctx, tau=None, q=None, m=None, k=None, qbar=None,
|
||||
g2=None, g3=None, omega1=None, omega2=None):
|
||||
r"""
|
||||
Evaluates the Klein j-invariant, which is a modular function defined for
|
||||
`\tau` in the upper half-plane as
|
||||
@@ -400,6 +467,11 @@ def kleinj(ctx, tau=None, **kwargs):
|
||||
An alternative, common notation is that of the j-function
|
||||
`j(\tau) = 1728 J(\tau)`.
|
||||
|
||||
The half-period ratio may be supplied directly as ``tau``, or through
|
||||
any of the argument forms accepted by :func:`~mpmath.qfrom`.
|
||||
The invariant can also be computed directly from Weierstrass invariants
|
||||
by providing both ``g2`` and ``g3``.
|
||||
|
||||
**Plots**
|
||||
|
||||
.. literalinclude :: /plots/kleinj.py
|
||||
@@ -465,7 +537,13 @@ def kleinj(ctx, tau=None, **kwargs):
|
||||
1264538.909475140509320227
|
||||
|
||||
"""
|
||||
q = ctx.qfrom(tau=tau, **kwargs)
|
||||
_validate_weierstrass_pair_args("kleinj", g2, g3, omega1, omega2)
|
||||
if g2 is not None:
|
||||
g2 = ctx.convert(g2)
|
||||
g3 = ctx.convert(g3)
|
||||
return g2**3 / (g2**3 - 27*g3**2)
|
||||
q = ctx.qfrom(q=q, m=m, k=k, tau=tau, qbar=qbar,
|
||||
omega1=omega1, omega2=omega2)
|
||||
t2 = ctx.jtheta(2,0,q)
|
||||
t3 = ctx.jtheta(3,0,q)
|
||||
t4 = ctx.jtheta(4,0,q)
|
||||
@@ -474,6 +552,79 @@ def kleinj(ctx, tau=None, **kwargs):
|
||||
return P/Q
|
||||
|
||||
|
||||
@defun_wrapped
|
||||
def kleinjinv(ctx, J):
|
||||
r"""
|
||||
Evaluates a branch of the inverse Klein j-invariant.
|
||||
|
||||
Given a value `J`, returns a half-period ratio `\tau` in the upper
|
||||
half-plane such that ``kleinj(tau)`` equals `J`, up to numerical error.
|
||||
Since ``kleinj`` is invariant under modular transformations, the inverse
|
||||
is multivalued; this function returns one representative.
|
||||
The branch is determined by principal square/cube roots in the formula
|
||||
below; the returned value is one modularly equivalent preimage and is not
|
||||
canonicalized to a fundamental domain.
|
||||
|
||||
The implementation uses the classical inverse construction via the
|
||||
modular lambda function, described as Method 1 in the Wikipedia article
|
||||
at https://en.wikipedia.org/wiki/J-invariant#Inverse_functions.
|
||||
|
||||
It first converts from the normalized invariant `J` to `j = 1728 J`,
|
||||
then solves
|
||||
|
||||
.. math ::
|
||||
|
||||
j = 256 (1 - x)^3 / x^2, \quad x = \lambda (1 - \lambda),
|
||||
|
||||
for `lambda`. The half-period ratio is then obtained from
|
||||
|
||||
.. math ::
|
||||
|
||||
\tau = i K(1 - \lambda) / K(\lambda),
|
||||
|
||||
where `K` is the complete elliptic integral. Using
|
||||
|
||||
.. math ::
|
||||
|
||||
K(m) = \pi / (2 \operatorname{AGM}(1, \sqrt{1 - m})),
|
||||
|
||||
this is evaluated by the arithmetic-geometric mean formula
|
||||
|
||||
.. math ::
|
||||
|
||||
\tau =
|
||||
i \operatorname{AGM}(1, \sqrt{1 - \lambda}) /
|
||||
\operatorname{AGM}(1, \sqrt{\lambda}).
|
||||
|
||||
Different root choices give different modularly equivalent branches.
|
||||
|
||||
**Examples**
|
||||
|
||||
>>> from mpmath import mp, kleinj, kleinjinv, chop
|
||||
>>> mp.dps = 25
|
||||
>>> mp.pretty = True
|
||||
>>> tau = 0.625+0.75j
|
||||
>>> chop(kleinj(kleinjinv(kleinj(tau))) - kleinj(tau))
|
||||
0.0
|
||||
>>> kleinjinv(1)
|
||||
(0.0 + 1.0j)
|
||||
|
||||
"""
|
||||
J = ctx.convert(J)
|
||||
if J == 0:
|
||||
return ctx.mpc(-ctx.mpf(1)/2, ctx.sqrt(3)/2)
|
||||
_j = 1728 * J
|
||||
sqrt_arg = 3*(1728*_j**2 - _j**3)
|
||||
exponent = ctx.mpf(1) / ctx.mpf(3)
|
||||
t = (-_j**3 + 2304*_j**2 - 884736*_j +
|
||||
12288*ctx.sqrt(sqrt_arg))**exponent
|
||||
x = ctx.mpf(1)/768*t + (1 - _j/768) - (1536*_j - _j**2) / (768*t)
|
||||
|
||||
lbd = (1 + ctx.sqrt(1 - 4*x)) / 2
|
||||
tau = ctx.j * ctx.agm(1, ctx.sqrt(1-lbd)) / ctx.agm(1, ctx.sqrt(lbd))
|
||||
return tau
|
||||
|
||||
|
||||
def RF_calc(ctx, x, y, z, r):
|
||||
if y == z: return RC_calc(ctx, x, y, r)
|
||||
if x == z: return RC_calc(ctx, y, x, r)
|
||||
@@ -1518,43 +1669,6 @@ def _eisenstein_G4_G6(ctx, tau):
|
||||
G6 = 2 * ctx.zeta(6) * E6
|
||||
return G4, G6
|
||||
|
||||
def _inverse_kleinj(ctx, J):
|
||||
"""
|
||||
Compute tau from Klein's J-invariant using the inverse j-function.
|
||||
See: https://en.wikipedia.org/wiki/J-invariant
|
||||
"""
|
||||
J = ctx.convert(J)
|
||||
_j = 1728 * J
|
||||
sqrt_arg = 3*(1728*_j**2 - _j**3)
|
||||
exponent = ctx.mpf(1) / ctx.mpf(3)
|
||||
t = (-_j**3 + 2304*_j**2 - 884736*_j +
|
||||
12288*ctx.sqrt(sqrt_arg))**exponent
|
||||
x = ctx.mpf(1)/768*t + (1 - _j/768) - (1536*_j - _j**2) / (768*t)
|
||||
|
||||
lbd = (1 + ctx.sqrt(1 - 4*x)) / 2
|
||||
tau = ctx.j * ctx.agm(1, ctx.sqrt(1-lbd)) / ctx.agm(1, ctx.sqrt(lbd))
|
||||
return tau
|
||||
|
||||
def _kleinj_from_g2g3(ctx, g2, g3):
|
||||
"""
|
||||
Klein's absolute invariant J from g2, g3.
|
||||
(Not the j one with 1728 factor)
|
||||
https://mathworld.wolfram.com/KleinsAbsoluteInvariant.html
|
||||
"""
|
||||
g2 = ctx.convert(g2)
|
||||
g3 = ctx.convert(g3)
|
||||
return 1 / (1 - 27*g3**2/g2**3)
|
||||
|
||||
def _tau_from_g(ctx, g2, g3):
|
||||
"""
|
||||
Compute tau (half-period ratio) from g2, g3.
|
||||
"""
|
||||
g2 = ctx.convert(g2)
|
||||
g3 = ctx.convert(g3)
|
||||
J = _kleinj_from_g2g3(ctx, g2, g3)
|
||||
tau = _inverse_kleinj(ctx, J)
|
||||
return tau
|
||||
|
||||
def _weierstrass_omega_tau(ctx, funcname, g2=None, g3=None, tau=None,
|
||||
omega1=None, omega2=None):
|
||||
"""
|
||||
@@ -1572,10 +1686,7 @@ def _weierstrass_omega_tau(ctx, funcname, g2=None, g3=None, tau=None,
|
||||
if omega1 is not None:
|
||||
omega1 = ctx.convert(omega1)
|
||||
omega2 = ctx.convert(omega2)
|
||||
tau = omega2 / omega1
|
||||
if ctx.im(tau) <= 0:
|
||||
raise ValueError("%s: omega ratio must be in upper half-plane" %
|
||||
funcname)
|
||||
tau = ctx.taufrom(omega1=omega1, omega2=omega2)
|
||||
return omega1, tau
|
||||
if tau is not None:
|
||||
tau = ctx.convert(tau)
|
||||
@@ -1636,7 +1747,7 @@ def weierhalfperiods(ctx, g2, g3):
|
||||
>>> chop(g2), chop(g3)
|
||||
(60.0, 140.0)
|
||||
>>> chop(omega2/omega1)
|
||||
(0.5 + 0.209032224450873j)
|
||||
(-0.5 + 0.209032224450873j)
|
||||
|
||||
"""
|
||||
with ctx.extraprec(10):
|
||||
@@ -1648,12 +1759,12 @@ def weierhalfperiods(ctx, g2, g3):
|
||||
ctx.gamma(ctx.mpf(1)/ctx.mpf(3))**3 / (4*ctx.pi))
|
||||
tau = ctx.mpc(ctx.mpf(1)/ctx.mpf(2), ctx.sqrt(3)/2)
|
||||
elif g3 == 0:
|
||||
tau = _tau_from_g(ctx, g2, g3)
|
||||
tau = ctx.taufrom(g2=g2, g3=g3)
|
||||
G4, G6 = _eisenstein_G4_G6(ctx, tau)
|
||||
omegaA = (ctx.j * (ctx.mpf(15)/(4*g2) * G4) **
|
||||
(ctx.mpf(1)/ctx.mpf(4)))
|
||||
else:
|
||||
tau = _tau_from_g(ctx, g2, g3)
|
||||
tau = ctx.taufrom(g2=g2, g3=g3)
|
||||
G4, G6 = _eisenstein_G4_G6(ctx, tau)
|
||||
omegaA = ctx.sqrt(g2/g3 * G6/G4 * ctx.mpf(7)/ctx.mpf(12))
|
||||
|
||||
@@ -1679,7 +1790,8 @@ def weierhalfperiods(ctx, g2, g3):
|
||||
|
||||
scale = max([ctx.one] + [abs(x) for x in [e1, e2, e3] + wps])
|
||||
tolerance = ctx.sqrt(ctx.eps) * scale
|
||||
if mae > tolerance: raise ValueError("weierhalfperiods: no convergence")
|
||||
if mae > tolerance:
|
||||
raise ValueError("weierhalfperiods: no convergence")
|
||||
|
||||
omega1, omega2 = [omegas[k] for k in index_combos[min_index]][:2]
|
||||
if ctx.im(omega2/omega1) <= 0:
|
||||
@@ -1935,7 +2047,7 @@ def weierpinv(ctx, p, g2=None, g3=None, tau=None, omega1=None, omega2=None,
|
||||
- `g2, g3`: elliptic invariants
|
||||
- `tau` or `omega1, omega2`: alternative parameterizations
|
||||
- `weierp_prime` (optional): derivative value used to choose the sign of
|
||||
the inverse
|
||||
the inverse
|
||||
|
||||
**Examples**
|
||||
|
||||
|
||||
@@ -14,12 +14,14 @@ Author of the first version: M.T. Taschuk
|
||||
import random
|
||||
|
||||
import pytest
|
||||
import mpmath.functions.elliptic as elliptic_functions
|
||||
|
||||
from mpmath import (cos, cosh, cot, coth, csc, csch, diff, ellipe, ellipfun,
|
||||
ellipk, ellippi, elliprc, elliprd, elliprf, elliprg,
|
||||
elliprj, eps, exp, gamma, inf, isnan, j, jtheta, kleinj,
|
||||
ldexp, ln2, mp, mpc, mpf, nan, nsum, pi, polyroots, qfrom,
|
||||
sec, sech, sin, sinh, sqrt, tan, tanh, weierhalfperiods,
|
||||
kleinjinv, kfrom, ldexp, ln2, mp, mpc, mpf, mfrom, nan,
|
||||
nsum, pi, polyroots, qbarfrom, qfrom, sec, sech, sin,
|
||||
sinh, sqrt, tan, tanh, taufrom, weierhalfperiods,
|
||||
weierinvariants, weierp, weierpinv, weierpprime,
|
||||
weiersigma, weierzeta)
|
||||
|
||||
@@ -772,6 +774,31 @@ def test_elliptic_integrals():
|
||||
def test_issue_238():
|
||||
assert isnan(qfrom(m=nan))
|
||||
|
||||
def test_argument_conversions_from_weierstrass_data():
|
||||
mp.dps = 30
|
||||
|
||||
half_periods = [
|
||||
(1, j/2),
|
||||
(1, mpf(1)/2 + 3*j/4),
|
||||
(mpf(3)/4 + j/4, -mpf(1)/5 + 9*j/10),
|
||||
]
|
||||
|
||||
for omega1, omega2 in half_periods:
|
||||
tau = omega2 / omega1
|
||||
g2, g3 = weierinvariants(omega1, omega2)
|
||||
tau_from_invariants = taufrom(g2=g2, g3=g3)
|
||||
|
||||
for func in [qfrom, qbarfrom, kfrom, mfrom]:
|
||||
assert mpc_ae(func(omega1=omega1, omega2=omega2),
|
||||
func(tau=tau), eps=eps*1000)
|
||||
assert mpc_ae(func(g2=g2, g3=g3),
|
||||
func(tau=tau_from_invariants), eps=eps*1000)
|
||||
|
||||
g2, g3 = weierinvariants(1, j/2)
|
||||
for func in [qfrom, qbarfrom, kfrom, mfrom]:
|
||||
pytest.raises(ValueError, lambda func=func: func(g2=g2))
|
||||
pytest.raises(ValueError, lambda func=func: func(omega1=1))
|
||||
|
||||
def test_issue_604():
|
||||
assert ellipe(pi, 1).ae('2.0')
|
||||
|
||||
@@ -865,6 +892,13 @@ def test_weierstrass_half_periods_high_precision():
|
||||
assert mpc_ae(g2_roundtrip, g2, eps=eps*10000)
|
||||
assert mpc_ae(g3_roundtrip, g3, eps=eps*10000)
|
||||
|
||||
def test_weierstrass_half_periods_no_convergence(monkeypatch):
|
||||
def bad_roots(ctx, omega1, omega2):
|
||||
return [ctx.mpf(10), ctx.mpf(20), ctx.mpf(30)]
|
||||
|
||||
monkeypatch.setattr(elliptic_functions, "_roots_from_omega", bad_roots)
|
||||
pytest.raises(ValueError, lambda: weierhalfperiods(0, 1))
|
||||
|
||||
def test_weierstrass_parameter_conversions_with_kleinj():
|
||||
mp.dps = 30
|
||||
|
||||
@@ -876,6 +910,64 @@ def test_weierstrass_parameter_conversions_with_kleinj():
|
||||
|
||||
assert mpc_ae(kleinj(tau), j_from_invariants, eps=eps*1000)
|
||||
assert mpc_ae(kleinj(recovered_tau), kleinj(tau), eps=eps*1000)
|
||||
assert mpc_ae(kleinj(taufrom(g2=g2, g3=g3)), kleinj(tau),
|
||||
eps=eps*1000)
|
||||
|
||||
def test_kleinj_from_weierstrass_invariants():
|
||||
mp.dps = 30
|
||||
|
||||
tau = 0.625 + 0.75j
|
||||
g2, g3 = weierinvariants(0.5, tau/2)
|
||||
|
||||
assert mpc_ae(kleinj(g2=g2, g3=g3), kleinj(tau), eps=eps*1000)
|
||||
assert kleinj(g2=0, g3=1).ae(0)
|
||||
assert kleinj(g2=1, g3=0).ae(1)
|
||||
pytest.raises(ValueError, lambda: kleinj(g2=g2))
|
||||
pytest.raises(ValueError, lambda: kleinj(g3=g3))
|
||||
|
||||
def test_kleinj_argument_conversions():
|
||||
mp.dps = 30
|
||||
|
||||
tau = 0.625 + 0.75j
|
||||
value = kleinj(tau)
|
||||
|
||||
assert mpc_ae(kleinj(q=qfrom(tau=tau)), value, eps=eps*1000)
|
||||
assert mpc_ae(kleinj(qbar=qbarfrom(tau=tau)), value, eps=eps*1000)
|
||||
assert mpc_ae(kleinj(k=kfrom(tau=tau)), value, eps=eps*1000)
|
||||
assert mpc_ae(kleinj(m=mfrom(tau=tau)), value, eps=eps*1000)
|
||||
assert mpc_ae(kleinj(omega1=1, omega2=tau), value, eps=eps*1000)
|
||||
pytest.raises(ValueError, lambda: kleinj(omega1=1))
|
||||
|
||||
def test_kleinjinv():
|
||||
mp.dps = 30
|
||||
|
||||
tau = 0.625 + 0.75j
|
||||
value = kleinj(tau)
|
||||
|
||||
assert mpc_ae(kleinj(kleinjinv(value)), value, eps=eps*1000)
|
||||
assert mpc_ae(kleinjinv(0), -0.5 + sqrt(3)*j/2,
|
||||
eps=eps*1000)
|
||||
assert mpc_ae(kleinjinv(1), j, eps=eps*1000)
|
||||
|
||||
def test_taufrom_weierstrass_invariants():
|
||||
mp.dps = 30
|
||||
|
||||
tau = 0.625 + 0.75j
|
||||
g2, g3 = weierinvariants(0.5, tau/2)
|
||||
recovered_tau = taufrom(g2=g2, g3=g3)
|
||||
|
||||
assert mpc_ae(kleinj(recovered_tau), kleinj(tau), eps=eps*1000)
|
||||
pytest.raises(ValueError, lambda: taufrom(g2=g2))
|
||||
pytest.raises(ValueError, lambda: taufrom(g3=g3))
|
||||
|
||||
def test_taufrom_half_periods():
|
||||
mp.dps = 30
|
||||
|
||||
assert mpc_ae(taufrom(omega1=1, omega2=0.5*j), 0.5*j,
|
||||
eps=eps*1000)
|
||||
pytest.raises(ValueError, lambda: taufrom(omega1=1))
|
||||
pytest.raises(ValueError, lambda: taufrom(omega2=j))
|
||||
pytest.raises(ValueError, lambda: taufrom(omega1=1, omega2=-j))
|
||||
|
||||
def test_weierstrass_half_period_values_are_cubic_roots():
|
||||
mp.dps = 30
|
||||
|
||||
Reference in New Issue
Block a user