diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,14 @@ # Changelog for `weierstrass-functions` +## 0.1.6.0 - 2023-10-19 + +- Increased the lower bound of the version of the 'jacobi-theta' dependency. + +- Now the Weierstrass functions are defined from the half-periods rather than +from the elliptic invariants. Each of them has a "prime" version defined from +the elliptic invariants: `weierstrassP` and `weierstrassP'`, etc. + + ## 0.1.5.0 - 2023-10-18 Increased the lower bound of the version of the 'jacobi-theta' dependency, diff --git a/src/Math/Eisenstein.hs b/src/Math/Eisenstein.hs --- a/src/Math/Eisenstein.hs +++ b/src/Math/Eisenstein.hs @@ -22,7 +22,7 @@ import Data.Complex ( Complex(..) ) import Internal ( (%^%) ) import Math.EllipticIntegrals ( ellipticF', ellipticE' ) -import Math.JacobiTheta ( jtheta2, jtheta3, jtheta4, jtheta1Dash0 ) +import Math.JacobiTheta ( jtheta2', jtheta3', jtheta4', jtheta1Dash0 ) i_ :: Complex Double @@ -34,9 +34,8 @@ -> Complex Double lambda tau = (j2 / j3) %^% 4 where - q = exp (i_ * pi * tau) - j2 = jtheta2 0 q - j3 = jtheta3 0 q + j2 = jtheta2' 0 tau + j3 = jtheta3' 0 tau @@ -67,9 +66,8 @@ eisensteinE2 tau = 6 / pi * ellE * j3 - j3 * j3 - j4 where - q = exp (i_ * pi * tau) - j3 = jtheta3 0 q %^% 2 - j4 = jtheta4 0 q %^% 4 + j3 = jtheta3' 0 tau %^% 2 + j4 = jtheta4' 0 tau %^% 4 ellE = ellipticE' 1e-14 (pi/2) (lambda tau) -- | Eisenstein series of weight 4 @@ -77,19 +75,15 @@ Complex Double -- ^ tau -> Complex Double eisensteinE4 tau = - (jtheta2 0 q %^% 8 + jtheta3 0 q %^% 8 + jtheta4 0 q %^% 8) / 2 - where - q = exp (i_ * pi * tau) + (jtheta2' 0 tau %^% 8 + jtheta3' 0 tau %^% 8 + jtheta4' 0 tau %^% 8) / 2 -- | Eisenstein series of weight 6 eisensteinE6 :: Complex Double -- ^ tau -> Complex Double eisensteinE6 tau = - (jtheta3 0 q %^% 12 + jtheta4 0 q %^% 12 - 3 * jtheta2 0 q %^% 8 - * (jtheta3 0 q %^% 4 + jtheta4 0 q %^% 4)) / 2 - where - q = exp (i_ * pi * tau) + (jtheta3' 0 tau %^% 12 + jtheta4' 0 tau %^% 12 - 3 * jtheta2' 0 tau %^% 8 + * (jtheta3' 0 tau %^% 4 + jtheta4' 0 tau %^% 4)) / 2 -- | Modular discriminant modularDiscriminant :: @@ -144,8 +138,7 @@ etaDedekind tau = exp (ipitau / 12) * j3 where ipitau = i_ * pi * tau - q = exp (3 * ipitau) - j3 = jtheta3 (pi / 2 * (tau + 1)) q + j3 = jtheta3' (pi / 2 * (tau + 1)) (3 * tau) -- | Third derivative at 0 of the first Jacobi theta function jtheta1DashDashDash0 :: diff --git a/src/Math/Weierstrass.hs b/src/Math/Weierstrass.hs --- a/src/Math/Weierstrass.hs +++ b/src/Math/Weierstrass.hs @@ -11,10 +11,15 @@ ( halfPeriods, ellipticInvariants, weierstrassP, + weierstrassP', weierstrassPdash, + weierstrassPdash', weierstrassPinv, + weierstrassPinv', weierstrassSigma, - weierstrassZeta + weierstrassSigma', + weierstrassZeta, + weierstrassZeta' ) where import Data.Complex ( Complex(..) ) import Internal ( (%^%) ) @@ -22,10 +27,10 @@ eisensteinE6, kleinJinv, jtheta1DashDashDash0 ) -import Math.JacobiTheta ( jtheta2, - jtheta3, - jtheta1, - jtheta4, +import Math.JacobiTheta ( jtheta1', + jtheta2', + jtheta3', + jtheta4', jtheta1Dash0, jtheta1Dash ) import Math.Gamma ( gamma ) @@ -81,9 +86,8 @@ Complex Double -> Complex Double -> (Complex Double, Complex Double) g_from_omega1_and_tau omega1 tau = (g2, g3) where - q = exp (i_ * pi * tau) - j2 = jtheta2 0 q - j3 = jtheta3 0 q + j2 = jtheta2' 0 tau + j3 = jtheta3' 0 tau j2pow4 = j2 %^% 4 j2pow8 = j2pow4 * j2pow4 j2pow12 = j2pow4 * j2pow8 @@ -106,91 +110,120 @@ weierstrassP_from_tau z tau = (pi * j2 * j3 * j4 / j1) %^% 2 - pi * pi * (j2 %^% 4 + j3 %^% 4) / 3 where - q = exp (i_ * pi * tau) - j2 = jtheta2 0 q - j3 = jtheta3 0 q + j2 = jtheta2' 0 tau + j3 = jtheta3' 0 tau z' = pi * z - j1 = jtheta1 z' q - j4 = jtheta4 z' q + j1 = jtheta1' z' tau + j4 = jtheta4' z' tau -weierstrassP_from_omega :: - Complex Double -> Complex Double -> Complex Double -> Complex Double -weierstrassP_from_omega z omega1 omega2 = +-- | Weierstrass p-function given the half-periods +weierstrassP :: + Complex Double -- ^ z + -> Complex Double -- ^ half-period omega1 + -> Complex Double -- ^ half-period omega2 + -> Complex Double +weierstrassP z omega1 omega2 = weierstrassP_from_tau (z / omega1 / 2) (omega2 / omega1) / (4 * omega1 * omega1) --- | Weierstrass p-function -weierstrassP :: +-- | Weierstrass p-function given the elliptic invariants +weierstrassP' :: Complex Double -- ^ z -> Complex Double -- ^ elliptic invariant g2 -> Complex Double -- ^ elliptic invariant g3 -> Complex Double -weierstrassP z g2 g3 = weierstrassP_from_omega z omega1 omega2 +weierstrassP' z g2 g3 = weierstrassP z omega1 omega2 where (omega1, omega2) = halfPeriods g2 g3 --- | Derivative of Weierstrass p-function +-- | Derivative of Weierstrass p-function given the half-periods weierstrassPdash :: Complex Double -- ^ z - -> Complex Double -- ^ elliptic invariant g2 - -> Complex Double -- ^ elliptic invariant g3 + -> Complex Double -- ^ half-period omega1 + -> Complex Double -- ^ half-period omega2 -> Complex Double -weierstrassPdash z g2 g3 = 2 / (w1 %^% 3) * j2 * j3 * j4 * f +weierstrassPdash z omega1 omega2 = 2 / (w1 %^% 3) * j2 * j3 * j4 * f where - (omega1, omega2) = halfPeriods g2 g3 w1 = 2 * omega1 / pi tau = omega2 / omega1 q = exp (i_ * pi * tau) z' = z / w1 - j1 = jtheta1 z' q - j2 = jtheta2 z' q - j3 = jtheta3 z' q - j4 = jtheta4 z' q + j1 = jtheta1' z' tau + j2 = jtheta2' z' tau + j3 = jtheta3' z' tau + j4 = jtheta4' z' tau j1dash = jtheta1Dash0 q - j2zero = jtheta2 0 q - j3zero = jtheta3 0 q - j4zero = jtheta4 0 q + j2zero = jtheta2' 0 tau + j3zero = jtheta3' 0 tau + j4zero = jtheta4' 0 tau f = j1dash %^% 3 / (j1 %^% 3 * j2zero * j3zero * j4zero) --- | Inverse of Weierstrass p-function -weierstrassPinv :: - Complex Double -- ^ w +-- | Derivative of Weierstrass p-function given the elliptic invariants +weierstrassPdash' :: + Complex Double -- ^ z -> Complex Double -- ^ elliptic invariant g2 -> Complex Double -- ^ elliptic invariant g3 -> Complex Double -weierstrassPinv w g2 g3 = carlsonRF' 1e-14 (w - e1) (w - e2) (w - e3) +weierstrassPdash' z g2 g3 = weierstrassPdash z omega1 omega2 where (omega1, omega2) = halfPeriods g2 g3 - e1 = weierstrassP omega1 g2 g3 - e2 = weierstrassP omega2 g2 g3 - e3 = weierstrassP (-omega1 - omega2) g2 g3 --- | Weierstrass sigma function -weierstrassSigma :: +-- | Inverse of Weierstrass p-function given the half-periods +weierstrassPinv :: + Complex Double -- ^ w + -> Complex Double -- ^ half-period omega1 + -> Complex Double -- ^ half-period omega2 + -> Complex Double +weierstrassPinv w omega1 omega2 = carlsonRF' 1e-14 (w - e1) (w - e2) (w - e3) + where + e1 = weierstrassP omega1 omega1 omega2 + e2 = weierstrassP omega2 omega1 omega2 + e3 = weierstrassP (-omega1 - omega2) omega1 omega2 + +-- | Inverse of Weierstrass p-function given the elliptic invariants +weierstrassPinv' :: Complex Double -- ^ z -> Complex Double -- ^ elliptic invariant g2 -> Complex Double -- ^ elliptic invariant g3 -> Complex Double -weierstrassSigma z g2 g3 = w1 * exp (h * z * z1 / pi) * j1 / j1dash +weierstrassPinv' z g2 g3 = weierstrassPinv z omega1 omega2 where (omega1, omega2) = halfPeriods g2 g3 + +-- | Weierstrass sigma function given the half-periods +weierstrassSigma :: + Complex Double -- ^ z + -> Complex Double -- ^ half-period omega1 + -> Complex Double -- ^ half-period omega2 + -> Complex Double +weierstrassSigma z omega1 omega2 = w1 * exp (h * z * z1 / pi) * j1 / j1dash + where tau = omega2 / omega1 q = exp (i_ * pi * tau) w1 = -2 * omega1 / pi z1 = z / w1 - j1 = jtheta1 z1 q + j1 = jtheta1' z1 tau j1dash = jtheta1Dash0 q h = - pi / (6 * w1) * jtheta1DashDashDash0 tau / j1dash --- | Weierstrass zeta function -weierstrassZeta :: +-- | Weierstrass sigma function given the elliptic invariants +weierstrassSigma' :: Complex Double -- ^ z -> Complex Double -- ^ elliptic invariant g2 -> Complex Double -- ^ elliptic invariant g3 -> Complex Double -weierstrassZeta z g2 g3 = - eta1 * z + p * lj1dash +weierstrassSigma' z g2 g3 = weierstrassSigma z omega1 omega2 where (omega1, omega2) = halfPeriods g2 g3 + +-- | Weierstrass zeta function given the half-periods +weierstrassZeta :: + Complex Double -- ^ z + -> Complex Double -- ^ half-period omega1 + -> Complex Double -- ^ half-period omega2 + -> Complex Double +weierstrassZeta z omega1 omega2 = - eta1 * z + p * lj1dash + where tau = omega2 / omega1 q = exp (i_ * pi * tau) w1 = - omega1 / pi @@ -198,4 +231,14 @@ j1dash = jtheta1Dash0 q eta1 = p * jtheta1DashDashDash0 tau / (6 * w1 * j1dash) pz = p * z - lj1dash = jtheta1Dash pz q / jtheta1 pz q + lj1dash = jtheta1Dash pz q / jtheta1' pz tau + +-- | Weierstrass zeta function given the elliptic invariants +weierstrassZeta' :: + Complex Double -- ^ z + -> Complex Double -- ^ elliptic invariant g2 + -> Complex Double -- ^ elliptic invariant g3 + -> Complex Double +weierstrassZeta' z g2 g3 = weierstrassZeta z omega1 omega2 + where + (omega1, omega2) = halfPeriods g2 g3 diff --git a/tests/Main.hs b/tests/Main.hs --- a/tests/Main.hs +++ b/tests/Main.hs @@ -14,11 +14,11 @@ import Test.Tasty.HUnit ( testCase ) import Math.Weierstrass ( halfPeriods, ellipticInvariants, - weierstrassP, - weierstrassPdash, - weierstrassPinv, - weierstrassSigma, - weierstrassZeta ) + weierstrassP', + weierstrassPdash', + weierstrassPinv', + weierstrassSigma', + weierstrassZeta' ) i_ :: Complex Double i_ = 0.0 :+ 1.0 @@ -97,14 +97,14 @@ let z = 0.1 :+ 0.1 g2 = 2 :+ 1 g3 = 2 :+ (-1) - obtained = weierstrassP z g2 g3 + obtained = weierstrassP' z g2 g3 expected = (-0.0010285443715) :+ (-49.9979857342848) assertApproxEqual "" 11 expected obtained, testCase "Equianharmonic case" $ do let omega2 = gamma (1/3) ** 3 / 4 / pi z0 = omega2 * (1 :+ (1 / sqrt 3)) - obtained = weierstrassP z0 0 1 + obtained = weierstrassP' z0 0 1 expected = 0 assertApproxEqual "" 13 obtained expected, @@ -112,8 +112,8 @@ let z = 1 :+ 1 g2 = 5 :+ 3 g3 = 2 :+ 7 - w = weierstrassP z g2 g3 - wdash = weierstrassPdash z g2 g3 + w = weierstrassP' z g2 g3 + wdash = weierstrassPdash' z g2 g3 left = wdash ** 2 right = 4 * w ** 3 - g2 * w - g3 assertApproxEqual "" 11 left right, @@ -122,8 +122,8 @@ let w = 0.1 :+ 1 g2 = 2 :+ 2 g3 = 0 :+ 3 - z = weierstrassPinv w g2 g3 - obtained = weierstrassP z g2 g3 + z = weierstrassPinv' w g2 g3 + obtained = weierstrassP' z g2 g3 expected = w assertApproxEqual "" 13 expected obtained, @@ -140,14 +140,14 @@ testCase "a value of weierstrassSigma" $ do let expected = 1.8646253716 :+ (-0.3066001355) - obtained = weierstrassSigma 2 1 (2 * i_) + obtained = weierstrassSigma' 2 1 (2 * i_) assertApproxEqual "" 10 expected obtained, testCase "a value of weierstrassZeta" $ do let g2 = 5 :+ 3 g3 = 5 :+ 3 expected = 0.802084165492408 :+ (-0.381791358666872) - obtained = weierstrassZeta (1 :+ 1) g2 g3 + obtained = weierstrassZeta' (1 :+ 1) g2 g3 assertApproxEqual "" 13 expected obtained ] diff --git a/weierstrass-functions.cabal b/weierstrass-functions.cabal --- a/weierstrass-functions.cabal +++ b/weierstrass-functions.cabal @@ -1,5 +1,5 @@ name: weierstrass-functions -version: 0.1.5.0 +version: 0.1.6.0 synopsis: Weierstrass Elliptic Functions description: Evaluation of Weierstrass elliptic functions and some related functions. homepage: https://github.com/stla/weierstrass-functions#readme @@ -20,7 +20,7 @@ , Math.Weierstrass other-modules: Internal build-depends: base >= 4.7 && < 5 - , jacobi-theta >= 0.2.1.1 + , jacobi-theta >= 0.2.2.0 , elliptic-integrals >= 0.1.0.0 , gamma >= 0.10.0.0 default-language: Haskell2010