weierstrass-functions 0.1.5.0 → 0.1.6.0

raw patch · 5 files changed

+122/−77 lines, 5 filesdep ~jacobi-thetaPVP ok

version bump matches the API change (PVP)

Dependency ranges changed: jacobi-theta

API changes (from Hackage documentation)

Files

• CHANGELOG.md +9/−0
• src/Math/Eisenstein.hs +9/−16
• src/Math/Weierstrass.hs +89/−46
• tests/Main.hs +13/−13
• weierstrass-functions.cabal +2/−2

@@ -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, 

@@ -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 :: 

@@ -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

@@ -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
 
   ]

@@ -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