weierstrass-functions 0.1.1.0 → 0.1.2.0
raw patch · 10 files changed
+596/−577 lines, 10 filessetup-changedPVP ok
version bump matches the API change (PVP)
API changes (from Hackage documentation)
Files
- CHANGELOG.md +15/−10
- LICENSE +30/−30
- README.md +7/−7
- Setup.hs +2/−2
- src/Internal.hs +5/−5
- src/Math/Eisenstein.hs +127/−122
- src/Math/Weierstrass.hs +191/−191
- tests/Approx.hs +15/−15
- tests/Main.hs +153/−144
- weierstrass-functions.cabal +51/−51
CHANGELOG.md view
@@ -1,10 +1,15 @@-# Changelog for `weierstrass-functions`--## 0.1.1.0 - 2023-03-02--Changes due to the update of 'jacobi-theta'.---## 0.1.0.0 - 2023-22-02--First release.+# Changelog for `weierstrass-functions` + +## 0.1.2.0 - 2023-08-24 + +Simpler and better implementation of `kleinJ`. + + +## 0.1.1.0 - 2023-03-02 + +Changes due to the update of 'jacobi-theta'. + + +## 0.1.0.0 - 2023-02-22 + +First release.
LICENSE view
@@ -1,30 +1,30 @@-Copyright Stéphane Laurent (c) 2023--All rights reserved.--Redistribution and use in source and binary forms, with or without-modification, are permitted provided that the following conditions are met:-- * Redistributions of source code must retain the above copyright- notice, this list of conditions and the following disclaimer.-- * Redistributions in binary form must reproduce the above- copyright notice, this list of conditions and the following- disclaimer in the documentation and/or other materials provided- with the distribution.-- * Neither the name of Stéphane Laurent nor the names of other- contributors may be used to endorse or promote products derived- from this software without specific prior written permission.--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS-"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT-LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR-A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT-OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,-SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT-LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,-DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY-THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT-(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE-OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.+Copyright Stéphane Laurent (c) 2023 + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Stéphane Laurent nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
README.md view
@@ -1,8 +1,8 @@-# weierstrass-functions--<!-- badges: start -->-[](https://github.com/stla/weierstrass-functions/actions/workflows/Stack-lts.yml)-[](https://github.com/stla/weierstrass-functions/actions/workflows/Stack-nightly.yml)-<!-- badges: end -->-+# weierstrass-functions + +<!-- badges: start --> +[](https://github.com/stla/weierstrass-functions/actions/workflows/Stack-lts.yml) +[](https://github.com/stla/weierstrass-functions/actions/workflows/Stack-nightly.yml) +<!-- badges: end --> + Evaluation of the Weierstrass elliptic functions and some related functions.
Setup.hs view
@@ -1,2 +1,2 @@-import Distribution.Simple-main = defaultMain+import Distribution.Simple +main = defaultMain
src/Internal.hs view
@@ -1,5 +1,5 @@-module Internal ((%^%)) where-import Data.Complex ( Complex (..) )--(%^%) :: Complex Double -> Int -> Complex Double-(%^%) z p = z ^^ p+module Internal ((%^%)) where +import Data.Complex ( Complex (..) ) + +(%^%) :: Complex Double -> Int -> Complex Double +(%^%) z p = z ^^ p
src/Math/Eisenstein.hs view
@@ -1,122 +1,127 @@-module Math.Eisenstein- ( lambda,- eisensteinE2,- eisensteinE4,- eisensteinE6,- kleinJ,- kleinJinv,- modularDiscriminant,- agm,- etaDedekind,- jtheta1DashDashDash0,- jtheta1Dash- ) where-import Data.Complex ( Complex(..) )-import Internal ( (%^%) )-import Math.EllipticIntegrals ( ellipticF', ellipticE' )-import Math.JacobiTheta ( jtheta2, jtheta3, jtheta4, jtheta1Dash )---i_ :: Complex Double-i_ = 0.0 :+ 1.0---- | Lambda modular function (square of elliptic modulus)-lambda :: - Complex Double -- ^ tau- -> Complex Double-lambda tau = (j2 / j3) %^% 4- where- q = exp (i_ * pi * tau)- j2 = jtheta2 0 q- j3 = jtheta3 0 q---- | Eisenstein series of weight 2-eisensteinE2 :: - Complex Double -- ^ tau- -> Complex Double-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- ellE = ellipticE' 1e-14 (pi/2) (lambda tau)---- | Eisenstein series of weight 4-eisensteinE4 :: - 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)---- | 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)---- | Modular discriminant-modularDiscriminant ::- Complex Double -- ^ tau- -> Complex Double-modularDiscriminant tau = - (eisensteinE4 tau %^% 3 - eisensteinE6 tau %^% 2) / 1728---- | Klein J-function-kleinJ :: - Complex Double -- ^ tau- -> Complex Double-kleinJ tau = - eisensteinE4 tau %^% 3 / modularDiscriminant tau---- | Arithmetic-geometric mean-agm :: - Complex Double -- ^ x - -> Complex Double -- ^ y- -> Complex Double-agm x y = - if x == 0 || y == 0 || x + y == 0- then 0- else pi/4 * (x + y) / ellipticF' 1e-15 (pi/2) (((x - y) / (x + y)) %^% 2)---- | Inverse Klein-J function-kleinJinv :: - Complex Double- -> Complex Double-kleinJinv j = - if j == 0- then 0.5 :+ (sqrt 3 / 2)- else i_ * agm 1 (sqrt(1 - lbd)) / agm 1 (sqrt lbd)- where- j2 = j * j- j3 = j2 * j- t = -j3 + 2304 * j2 + 12288 * sqrt(3 * (1728 * j2 - j3)) - 884736 * j- u = t ** (1/3)- x = 4 + (u - j) / 192 - (1536 * j - j2) / (192 * u)- lbd = -(-1 - sqrt(1 - x)) / 2---- | Dedekind eta function-etaDedekind ::- Complex Double -- ^ tau- -> Complex Double-etaDedekind tau = exp (ipitau / 12) * j3- where- ipitau = i_ * pi * tau- q = exp (3 * ipitau)- j3 = jtheta3 (pi / 2 * (tau + 1)) q---- | Third derivative at 0 of the first Jacobi theta function-jtheta1DashDashDash0 :: - Complex Double -- ^ tau- -> Complex Double-jtheta1DashDashDash0 tau = - jtheta1Dash 0 q * eisensteinE2 tau - where- q = exp (i_ * pi * tau)--+module Math.Eisenstein + ( lambda, + eisensteinE2, + eisensteinE4, + eisensteinE6, + kleinJ, + kleinJinv, + modularDiscriminant, + agm, + etaDedekind, + jtheta1DashDashDash0, + jtheta1Dash + ) where +import Data.Complex ( Complex(..) ) +import Internal ( (%^%) ) +import Math.EllipticIntegrals ( ellipticF', ellipticE' ) +import Math.JacobiTheta ( jtheta2, jtheta3, jtheta4, jtheta1Dash ) + + +i_ :: Complex Double +i_ = 0.0 :+ 1.0 + +-- | Lambda modular function (square of elliptic modulus) +lambda :: + Complex Double -- ^ tau + -> Complex Double +lambda tau = (j2 / j3) %^% 4 + where + q = exp (i_ * pi * tau) + j2 = jtheta2 0 q + j3 = jtheta3 0 q + +-- | Eisenstein series of weight 2 +eisensteinE2 :: + Complex Double -- ^ tau + -> Complex Double +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 + ellE = ellipticE' 1e-14 (pi/2) (lambda tau) + +-- | Eisenstein series of weight 4 +eisensteinE4 :: + 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) + +-- | 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) + +-- | Modular discriminant +modularDiscriminant :: + Complex Double -- ^ tau + -> Complex Double +modularDiscriminant tau = + (eisensteinE4 tau %^% 3 - eisensteinE6 tau %^% 2) / 1728 + +-- | Klein J-function +kleinJ :: + Complex Double -- ^ tau + -> Complex Double +kleinJ tau = + let + lbd = lambda(tau) + x = lbd * (1 - lbd) + in + 256 * (1 - x) %^% 3 / x %^% 2 + --eisensteinE4 tau %^% 3 / modularDiscriminant tau + +-- | Arithmetic-geometric mean +agm :: + Complex Double -- ^ x + -> Complex Double -- ^ y + -> Complex Double +agm x y = + if x == 0 || y == 0 || x + y == 0 + then 0 + else pi/4 * (x + y) / ellipticF' 1e-15 (pi/2) (((x - y) / (x + y)) %^% 2) + +-- | Inverse Klein-J function +kleinJinv :: + Complex Double + -> Complex Double +kleinJinv j = + if j == 0 + then 0.5 :+ (sqrt 3 / 2) + else i_ * agm 1 (sqrt(1 - lbd)) / agm 1 (sqrt lbd) + where + j2 = j * j + j3 = j2 * j + t = -j3 + 2304 * j2 + 12288 * sqrt(3 * (1728 * j2 - j3)) - 884736 * j + u = t ** (1/3) + x = 4 + (u - j) / 192 - (1536 * j - j2) / (192 * u) + lbd = -(-1 - sqrt(1 - x)) / 2 + +-- | Dedekind eta function +etaDedekind :: + Complex Double -- ^ tau + -> Complex Double +etaDedekind tau = exp (ipitau / 12) * j3 + where + ipitau = i_ * pi * tau + q = exp (3 * ipitau) + j3 = jtheta3 (pi / 2 * (tau + 1)) q + +-- | Third derivative at 0 of the first Jacobi theta function +jtheta1DashDashDash0 :: + Complex Double -- ^ tau + -> Complex Double +jtheta1DashDashDash0 tau = - jtheta1Dash 0 q * eisensteinE2 tau + where + q = exp (i_ * pi * tau) + +
src/Math/Weierstrass.hs view
@@ -1,191 +1,191 @@-module Math.Weierstrass- ( halfPeriods,- ellipticInvariants,- weierstrassP,- weierstrassPdash,- weierstrassPinv,- weierstrassSigma,- weierstrassZeta- ) where-import Data.Complex ( Complex(..) )-import Internal ( (%^%) )-import Math.Eisenstein ( eisensteinE4, - eisensteinE6, - kleinJinv, - jtheta1DashDashDash0 ) -import Math.JacobiTheta ( jtheta2, - jtheta3, - jtheta1, - jtheta4,- jtheta1Dash )-import Math.Gamma ( gamma )-import Math.EllipticIntegrals ( carlsonRF' )----i_ :: Complex Double-i_ = 0.0 :+ 1.0--eisensteinG4 :: Complex Double -> Complex Double-eisensteinG4 tau = pi %^% 4 / 45 * eisensteinE4 tau--eisensteinG6_over_eisensteinG4 :: Complex Double -> Complex Double-eisensteinG6_over_eisensteinG4 tau = - 2 * pi * pi / 21 * eisensteinE6 tau / eisensteinE4 tau--omega1_and_tau :: - Complex Double -> Complex Double -> (Complex Double, Complex Double)-omega1_and_tau g2 g3 = (omega1, tau)- where- (omega1, tau) - | g2 == 0 = - (- gamma (1/3) %^% 3 / (4 * pi * g3 ** (1/6)),- 0.5 :+ (sqrt 3 / 2)- )- | g3 == 0 = - (- i_ * sqrt(sqrt(3.75 * eisensteinG4 tau' / g2)),- tau'- )- | otherwise = - (- sqrt(7 * g2 * eisensteinG6_over_eisensteinG4 tau' / (12 * g3)),- tau'- )- where - g2cube = g2 * g2 * g2- j = 1728 * g2cube / (g2cube - 27 * g3 * g3)- tau' = kleinJinv j---- | Half-periods from elliptic invariants.-halfPeriods :: - Complex Double -- ^ g2- -> Complex Double -- ^ g3- -> (Complex Double, Complex Double) -- ^ omega1, omega2-halfPeriods g2 g3 = (omega1, tau * omega1)- where- (omega1, tau) = omega1_and_tau g2 g3--g_from_omega1_and_tau :: - 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- j2pow4 = j2 %^% 4- j2pow8 = j2pow4 * j2pow4- j2pow12 = j2pow4 * j2pow8- j3pow4 = j3 %^% 4- j3pow8 = j3pow4 * j3pow4- j3pow12 = j3pow4 * j3pow8- g2 = 4/3 * (pi / 2 / omega1) %^% 4 * (j2pow8 - j2pow4 * j3pow4 + j3pow8)- g3 = 8/27 * (pi / 2 / omega1) %^% 6 *- (j2pow12 - ((1.5 * j2pow8 * j3pow4) + (1.5 * j2pow4 * j3pow8)) + j3pow12)---- | Elliptic invariants from half-periods.-ellipticInvariants :: - Complex Double -- ^ omega1- -> Complex Double -- ^ omega2- -> (Complex Double, Complex Double) -- ^ g2, g3-ellipticInvariants omega1 omega2 = - g_from_omega1_and_tau omega1 (omega2 / omega1)--weierstrassP_from_tau :: Complex Double -> Complex Double -> Complex Double-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- z' = pi * z- j1 = jtheta1 z' q- j4 = jtheta4 z' q--weierstrassP_from_omega :: - Complex Double -> Complex Double -> Complex Double -> Complex Double-weierstrassP_from_omega z omega1 omega2 = - weierstrassP_from_tau - (z / omega1 / 2) (omega2 / omega1) / (4 * omega1 * omega1)---- | Weierstrass p-function-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- where- (omega1, omega2) = halfPeriods g2 g3---- | Derivative of Weierstrass p-function-weierstrassPdash ::- Complex Double -- ^ z- -> Complex Double -- ^ elliptic invariant g2- -> Complex Double -- ^ elliptic invariant g3- -> Complex Double-weierstrassPdash z g2 g3 = 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- j1dash = jtheta1Dash 0 q- j2zero = jtheta2 0 q- j3zero = jtheta3 0 q- j4zero = jtheta4 0 q- f = j1dash %^% 3 / (j1 %^% 3 * j2zero * j3zero * j4zero)---- | Inverse of Weierstrass p-function-weierstrassPinv ::- Complex Double -- ^ w- -> 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)- 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 ::- 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- where- (omega1, omega2) = halfPeriods g2 g3- tau = omega2 / omega1- q = exp (i_ * pi * tau)- w1 = -2 * omega1 / pi- z1 = z / w1- j1 = jtheta1 z1 q- j1dash = jtheta1Dash 0 q- h = - pi / (6 * w1) * jtheta1DashDashDash0 tau / j1dash---- | Weierstrass zeta function-weierstrassZeta ::- Complex Double -- ^ z- -> Complex Double -- ^ elliptic invariant g2- -> Complex Double -- ^ elliptic invariant g3- -> Complex Double-weierstrassZeta z g2 g3 = - eta1 * z + p * lj1dash- where- (omega1, omega2) = halfPeriods g2 g3- tau = omega2 / omega1- q = exp (i_ * pi * tau)- w1 = - omega1 / pi- p = 0.5 / w1- j1dash = jtheta1Dash 0 q- eta1 = p * jtheta1DashDashDash0 tau / (6 * w1 * j1dash)- pz = p * z- lj1dash = jtheta1Dash pz q / jtheta1 pz q+module Math.Weierstrass + ( halfPeriods, + ellipticInvariants, + weierstrassP, + weierstrassPdash, + weierstrassPinv, + weierstrassSigma, + weierstrassZeta + ) where +import Data.Complex ( Complex(..) ) +import Internal ( (%^%) ) +import Math.Eisenstein ( eisensteinE4, + eisensteinE6, + kleinJinv, + jtheta1DashDashDash0 ) +import Math.JacobiTheta ( jtheta2, + jtheta3, + jtheta1, + jtheta4, + jtheta1Dash ) +import Math.Gamma ( gamma ) +import Math.EllipticIntegrals ( carlsonRF' ) + + + +i_ :: Complex Double +i_ = 0.0 :+ 1.0 + +eisensteinG4 :: Complex Double -> Complex Double +eisensteinG4 tau = pi %^% 4 / 45 * eisensteinE4 tau + +eisensteinG6_over_eisensteinG4 :: Complex Double -> Complex Double +eisensteinG6_over_eisensteinG4 tau = + 2 * pi * pi / 21 * eisensteinE6 tau / eisensteinE4 tau + +omega1_and_tau :: + Complex Double -> Complex Double -> (Complex Double, Complex Double) +omega1_and_tau g2 g3 = (omega1, tau) + where + (omega1, tau) + | g2 == 0 = + ( + gamma (1/3) %^% 3 / (4 * pi * g3 ** (1/6)), + 0.5 :+ (sqrt 3 / 2) + ) + | g3 == 0 = + ( + i_ * sqrt(sqrt(3.75 * eisensteinG4 tau' / g2)), + tau' + ) + | otherwise = + ( + sqrt(7 * g2 * eisensteinG6_over_eisensteinG4 tau' / (12 * g3)), + tau' + ) + where + g2cube = g2 * g2 * g2 + j = 1728 * g2cube / (g2cube - 27 * g3 * g3) + tau' = kleinJinv j + +-- | Half-periods from elliptic invariants. +halfPeriods :: + Complex Double -- ^ g2 + -> Complex Double -- ^ g3 + -> (Complex Double, Complex Double) -- ^ omega1, omega2 +halfPeriods g2 g3 = (omega1, tau * omega1) + where + (omega1, tau) = omega1_and_tau g2 g3 + +g_from_omega1_and_tau :: + 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 + j2pow4 = j2 %^% 4 + j2pow8 = j2pow4 * j2pow4 + j2pow12 = j2pow4 * j2pow8 + j3pow4 = j3 %^% 4 + j3pow8 = j3pow4 * j3pow4 + j3pow12 = j3pow4 * j3pow8 + g2 = 4/3 * (pi / 2 / omega1) %^% 4 * (j2pow8 - j2pow4 * j3pow4 + j3pow8) + g3 = 8/27 * (pi / 2 / omega1) %^% 6 * + (j2pow12 - ((1.5 * j2pow8 * j3pow4) + (1.5 * j2pow4 * j3pow8)) + j3pow12) + +-- | Elliptic invariants from half-periods. +ellipticInvariants :: + Complex Double -- ^ omega1 + -> Complex Double -- ^ omega2 + -> (Complex Double, Complex Double) -- ^ g2, g3 +ellipticInvariants omega1 omega2 = + g_from_omega1_and_tau omega1 (omega2 / omega1) + +weierstrassP_from_tau :: Complex Double -> Complex Double -> Complex Double +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 + z' = pi * z + j1 = jtheta1 z' q + j4 = jtheta4 z' q + +weierstrassP_from_omega :: + Complex Double -> Complex Double -> Complex Double -> Complex Double +weierstrassP_from_omega z omega1 omega2 = + weierstrassP_from_tau + (z / omega1 / 2) (omega2 / omega1) / (4 * omega1 * omega1) + +-- | Weierstrass p-function +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 + where + (omega1, omega2) = halfPeriods g2 g3 + +-- | Derivative of Weierstrass p-function +weierstrassPdash :: + Complex Double -- ^ z + -> Complex Double -- ^ elliptic invariant g2 + -> Complex Double -- ^ elliptic invariant g3 + -> Complex Double +weierstrassPdash z g2 g3 = 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 + j1dash = jtheta1Dash 0 q + j2zero = jtheta2 0 q + j3zero = jtheta3 0 q + j4zero = jtheta4 0 q + f = j1dash %^% 3 / (j1 %^% 3 * j2zero * j3zero * j4zero) + +-- | Inverse of Weierstrass p-function +weierstrassPinv :: + Complex Double -- ^ w + -> 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) + 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 :: + 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 + where + (omega1, omega2) = halfPeriods g2 g3 + tau = omega2 / omega1 + q = exp (i_ * pi * tau) + w1 = -2 * omega1 / pi + z1 = z / w1 + j1 = jtheta1 z1 q + j1dash = jtheta1Dash 0 q + h = - pi / (6 * w1) * jtheta1DashDashDash0 tau / j1dash + +-- | Weierstrass zeta function +weierstrassZeta :: + Complex Double -- ^ z + -> Complex Double -- ^ elliptic invariant g2 + -> Complex Double -- ^ elliptic invariant g3 + -> Complex Double +weierstrassZeta z g2 g3 = - eta1 * z + p * lj1dash + where + (omega1, omega2) = halfPeriods g2 g3 + tau = omega2 / omega1 + q = exp (i_ * pi * tau) + w1 = - omega1 / pi + p = 0.5 / w1 + j1dash = jtheta1Dash 0 q + eta1 = p * jtheta1DashDashDash0 tau / (6 * w1 * j1dash) + pz = p * z + lj1dash = jtheta1Dash pz q / jtheta1 pz q
tests/Approx.hs view
@@ -1,15 +1,15 @@-module Approx (assertApproxEqual) where-import Data.Complex ( imagPart, realPart, Complex(..) )-import Test.Tasty.HUnit ( Assertion, assertEqual )---- round x to n digits-approx0 :: Int -> Double -> Double-approx0 n x = fromInteger (round $ x * (10^n)) / (10.0^^n)---- round z to n digits-approx :: Int -> Complex Double -> Complex Double-approx n z = approx0 n (realPart z) :+ approx0 n (imagPart z)--assertApproxEqual :: String -> Int -> Complex Double -> Complex Double -> Assertion-assertApproxEqual prefix n z1 z2 = - assertEqual prefix (approx n z1) (approx n z2)+module Approx (assertApproxEqual) where +import Data.Complex ( imagPart, realPart, Complex(..) ) +import Test.Tasty.HUnit ( Assertion, assertEqual ) + +-- round x to n digits +approx0 :: Int -> Double -> Double +approx0 n x = fromInteger (round $ x * (10^n)) / (10.0^^n) + +-- round z to n digits +approx :: Int -> Complex Double -> Complex Double +approx n z = approx0 n (realPart z) :+ approx0 n (imagPart z) + +assertApproxEqual :: String -> Int -> Complex Double -> Complex Double -> Assertion +assertApproxEqual prefix n z1 z2 = + assertEqual prefix (approx n z1) (approx n z2)
tests/Main.hs view
@@ -1,144 +1,153 @@-module Main where-import Approx ( assertApproxEqual )-import Data.Complex ( Complex(..) )-import Math.Eisenstein ( eisensteinE4,- eisensteinE6,- kleinJ,- agm,- kleinJinv, - etaDedekind,- lambda )-import Math.Gamma ( gamma )-import Test.Tasty ( defaultMain, testGroup )-import Test.Tasty.HUnit ( testCase )-import Math.Weierstrass ( halfPeriods, - ellipticInvariants,- weierstrassP,- weierstrassPdash,- weierstrassPinv,- weierstrassSigma,- weierstrassZeta )--i_ :: Complex Double-i_ = 0.0 :+ 1.0--tau1 :: Complex Double -tau1 = i_--tau2 :: Complex Double -tau2 = i_ / 10.0--tau3 :: Complex Double -tau3 = 2.0 :+ 2.0--main :: IO ()-main = defaultMain $- testGroup "Tests"- [ - testCase "E4 is modular - condition 1" $ do- let e4_tau = eisensteinE4 tau1 - e4_taup1 = eisensteinE4 (tau1 + 1)- assertApproxEqual "" 12 e4_tau e4_taup1,-- testCase "E4 is modular - condition 2" $ do- let e4 = eisensteinE4 (-1 / tau2) - e4' = tau2**4 * eisensteinE4 tau2- assertApproxEqual "" 12 e4 e4',-- testCase "E6 is modular - condition 1" $ do- let e6_tau = eisensteinE6 tau2 - e6_taup1 = eisensteinE6 (tau2 + 1)- assertApproxEqual "" 7 e6_tau e6_taup1,-- testCase "E6 is modular - condition 2" $ do- let e6 = eisensteinE6 (-1 / tau3) - e6' = tau3**6 * eisensteinE6 tau3- assertApproxEqual "" 10 e6 e6',-- testCase "a value of Klein J-function" $ do- let expected = 66**3- obtained = kleinJ (2 * i_)- assertApproxEqual "" 7 expected obtained,-- testCase "a value of agm" $ do- let expected = 2 * pi ** 1.5 * sqrt 2 / gamma 0.25 ** 2- obtained = agm 1 (sqrt 2)- assertApproxEqual "" 14 expected obtained,-- testCase "kleinJ o kleinJinv = id" $ do- let expected = 0.2 :+ 0.2- obtained = kleinJ (kleinJinv (0.2 :+ 0.2))- assertApproxEqual "" 12 expected obtained,-- testCase "Elliptic invariants - 1/2" $ do- let g2 = (-7) :+ 9- g3 = 5 :+ 3- (omega1, omega2) = halfPeriods g2 g3- (g2', _) = ellipticInvariants omega1 omega2- assertApproxEqual "" 12 g2 g2',-- testCase "Elliptic invariants - 2/2" $ do- let g2 = (-7) :+ 9- g3 = 5 :+ 3- (omega1, omega2) = halfPeriods g2 g3- (_, g3') = ellipticInvariants omega1 omega2- assertApproxEqual "" 12 g3 g3',-- testCase "a value of weierstrassP" $ do- let z = 0.1 :+ 0.1- g2 = 2 :+ 1- g3 = 2 :+ (-1)- 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- expected = 0- assertApproxEqual "" 13 obtained expected,-- testCase "Differential equation" $ do- let z = 1 :+ 1- g2 = 2 :+ 1- g3 = 2 :+ (-1)- w = weierstrassP z g2 g3- wdash = weierstrassPdash z g2 g3- left = wdash ** 2- right = 4 * w ** 3 - g2 * w - g3- assertApproxEqual "" 11 left right,-- testCase "weierstrassPinv works" $ do- let w = 0.1 :+ 1- g2 = 2 :+ 2- g3 = 0 :+ 3- z = weierstrassPinv w g2 g3- obtained = weierstrassP z g2 g3- expected = w- assertApproxEqual "" 13 expected obtained,-- testCase "a value of Dedekind eta" $ do- let expected = gamma 0.25 / 2 ** (11/8) / pi ** 0.75- obtained = etaDedekind (2 * i_)- assertApproxEqual "" 14 expected obtained,-- testCase "lambda modular identity" $ do- let x = sqrt 2- expected = 1- obtained = lambda (i_ * x) + lambda (i_ / x)- assertApproxEqual "" 14 expected obtained,-- testCase "a value of weierstrassSigma" $ do- let expected = 1.8646253716 :+ (-0.3066001355)- 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- assertApproxEqual "" 13 expected obtained-- ]+module Main where +import Approx ( assertApproxEqual ) +import Data.Complex ( Complex(..) ) +import Math.Eisenstein ( eisensteinE4, + eisensteinE6, + modularDiscriminant, + kleinJ, + agm, + kleinJinv, + etaDedekind, + lambda ) +import Math.Gamma ( gamma ) +import Test.Tasty ( defaultMain, testGroup ) +import Test.Tasty.HUnit ( testCase ) +import Math.Weierstrass ( halfPeriods, + ellipticInvariants, + weierstrassP, + weierstrassPdash, + weierstrassPinv, + weierstrassSigma, + weierstrassZeta ) + +i_ :: Complex Double +i_ = 0.0 :+ 1.0 + +tau1 :: Complex Double +tau1 = i_ + +tau2 :: Complex Double +tau2 = i_ / 10.0 + +tau3 :: Complex Double +tau3 = 2.0 :+ 2.0 + +tau4 :: Complex Double +tau4 = 0.2 :+ 0.2 + +main :: IO () +main = defaultMain $ + testGroup "Tests" + [ + testCase "E4 is modular - condition 1" $ do + let e4_tau = eisensteinE4 tau1 + e4_taup1 = eisensteinE4 (tau1 + 1) + assertApproxEqual "" 12 e4_tau e4_taup1, + + testCase "E4 is modular - condition 2" $ do + let e4 = eisensteinE4 (-1 / tau2) + e4' = tau2**4 * eisensteinE4 tau2 + assertApproxEqual "" 12 e4 e4', + + testCase "E6 is modular - condition 1" $ do + let e6_tau = eisensteinE6 tau2 + e6_taup1 = eisensteinE6 (tau2 + 1) + assertApproxEqual "" 7 e6_tau e6_taup1, + + testCase "E6 is modular - condition 2" $ do + let e6 = eisensteinE6 (-1 / tau3) + e6' = tau3**6 * eisensteinE6 tau3 + assertApproxEqual "" 10 e6 e6', + + testCase "a value of Klein J-function" $ do + let expected = 66**3 + obtained = kleinJ (2 * i_) + assertApproxEqual "" 7 expected obtained, + + testCase "a value of agm" $ do + let expected = 2 * pi ** 1.5 * sqrt 2 / gamma 0.25 ** 2 + obtained = agm 1 (sqrt 2) + assertApproxEqual "" 14 expected obtained, + + testCase "kleinJ o kleinJinv = id" $ do + let expected = 0.2 :+ 0.2 + obtained = kleinJ (kleinJinv (0.2 :+ 0.2)) + assertApproxEqual "" 12 expected obtained, + + testCase "kleinJ - alternative expression" $ do + let k = kleinJ tau3 + k' = (eisensteinE4 tau3)**3 / modularDiscriminant tau3 + assertApproxEqual "" 4 k k', + + testCase "Elliptic invariants - 1/2" $ do + let g2 = (-7) :+ 9 + g3 = 5 :+ 3 + (omega1, omega2) = halfPeriods g2 g3 + (g2', _) = ellipticInvariants omega1 omega2 + assertApproxEqual "" 12 g2 g2', + + testCase "Elliptic invariants - 2/2" $ do + let g2 = (-7) :+ 9 + g3 = 5 :+ 3 + (omega1, omega2) = halfPeriods g2 g3 + (_, g3') = ellipticInvariants omega1 omega2 + assertApproxEqual "" 12 g3 g3', + + testCase "a value of weierstrassP" $ do + let z = 0.1 :+ 0.1 + g2 = 2 :+ 1 + g3 = 2 :+ (-1) + 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 + expected = 0 + assertApproxEqual "" 13 obtained expected, + + testCase "Differential equation" $ do + let z = 1 :+ 1 + g2 = 2 :+ 1 + g3 = 2 :+ (-1) + w = weierstrassP z g2 g3 + wdash = weierstrassPdash z g2 g3 + left = wdash ** 2 + right = 4 * w ** 3 - g2 * w - g3 + assertApproxEqual "" 11 left right, + + testCase "weierstrassPinv works" $ do + let w = 0.1 :+ 1 + g2 = 2 :+ 2 + g3 = 0 :+ 3 + z = weierstrassPinv w g2 g3 + obtained = weierstrassP z g2 g3 + expected = w + assertApproxEqual "" 13 expected obtained, + + testCase "a value of Dedekind eta" $ do + let expected = gamma 0.25 / 2 ** (11/8) / pi ** 0.75 + obtained = etaDedekind (2 * i_) + assertApproxEqual "" 14 expected obtained, + + testCase "lambda modular identity" $ do + let x = sqrt 2 + expected = 1 + obtained = lambda (i_ * x) + lambda (i_ / x) + assertApproxEqual "" 14 expected obtained, + + testCase "a value of weierstrassSigma" $ do + let expected = 1.8646253716 :+ (-0.3066001355) + 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 + assertApproxEqual "" 13 expected obtained + + ]
weierstrass-functions.cabal view
@@ -1,51 +1,51 @@-name: weierstrass-functions-version: 0.1.1.0-synopsis: Weierstrass Elliptic Functions-description: Evaluation of Weierstrass elliptic functions and some related functions.-homepage: https://github.com/stla/weierstrass-functions#readme-license: BSD3-license-file: LICENSE-author: Stéphane Laurent-maintainer: laurent_step@outlook.fr-copyright: 2023 Stéphane Laurent-category: Math, Numeric-build-type: Simple-extra-source-files: README.md- CHANGELOG.md-cabal-version: >=1.10--library- hs-source-dirs: src- exposed-modules: Math.Eisenstein- , Math.Weierstrass- other-modules: Internal- build-depends: base >= 4.7 && < 5- , jacobi-theta >= 0.1.2.0- , elliptic-integrals >= 0.1.0.0- , gamma >= 0.10.0.0- default-language: Haskell2010- ghc-options: -Wall- -Wcompat- -Widentities- -Wincomplete-record-updates- -Wincomplete-uni-patterns- -Wmissing-export-lists- -Wmissing-home-modules- -Wpartial-fields- -Wredundant-constraints--test-suite unit-tests- type: exitcode-stdio-1.0- main-is: Main.hs- hs-source-dirs: tests/- other-modules: Approx- Build-Depends: base >= 4.7 && < 5- , tasty- , tasty-hunit- , weierstrass-functions- , gamma >= 0.10.0.0- Default-Language: Haskell2010--source-repository head- type: git- location: https://github.com/stla/weierstrass-functions+name: weierstrass-functions +version: 0.1.2.0 +synopsis: Weierstrass Elliptic Functions +description: Evaluation of Weierstrass elliptic functions and some related functions. +homepage: https://github.com/stla/weierstrass-functions#readme +license: BSD3 +license-file: LICENSE +author: Stéphane Laurent +maintainer: laurent_step@outlook.fr +copyright: 2023 Stéphane Laurent +category: Math, Numeric +build-type: Simple +extra-source-files: README.md + CHANGELOG.md +cabal-version: >=1.10 + +library + hs-source-dirs: src + exposed-modules: Math.Eisenstein + , Math.Weierstrass + other-modules: Internal + build-depends: base >= 4.7 && < 5 + , jacobi-theta >= 0.1.2.0 + , elliptic-integrals >= 0.1.0.0 + , gamma >= 0.10.0.0 + default-language: Haskell2010 + ghc-options: -Wall + -Wcompat + -Widentities + -Wincomplete-record-updates + -Wincomplete-uni-patterns + -Wmissing-export-lists + -Wmissing-home-modules + -Wpartial-fields + -Wredundant-constraints + +test-suite unit-tests + type: exitcode-stdio-1.0 + main-is: Main.hs + hs-source-dirs: tests/ + other-modules: Approx + Build-Depends: base >= 4.7 && < 5 + , tasty + , tasty-hunit + , weierstrass-functions + , gamma >= 0.10.0.0 + Default-Language: Haskell2010 + +source-repository head + type: git + location: https://github.com/stla/weierstrass-functions