weierstrass-functions (empty) → 0.1.0.0
raw patch · 10 files changed
+571/−0 lines, 10 filesdep +basedep +elliptic-integralsdep +gammasetup-changed
Dependencies added: base, elliptic-integrals, gamma, jacobi-theta, tasty, tasty-hunit, weierstrass-functions
Files
- CHANGELOG.md +6/−0
- LICENSE +30/−0
- README.md +8/−0
- Setup.hs +2/−0
- src/Internal.hs +5/−0
- src/Math/Eisenstein.hs +119/−0
- src/Math/Weierstrass.hs +191/−0
- tests/Approx.hs +15/−0
- tests/Main.hs +144/−0
- weierstrass-functions.cabal +51/−0
+ CHANGELOG.md view
@@ -0,0 +1,6 @@+# Changelog for `weierstrass-functions`+++## 0.1.0.0 - 2023-02-22++First release.
+ LICENSE view
@@ -0,0 +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.
+ README.md view
@@ -0,0 +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 -->++Evaluation of the Weierstrass elliptic functions and some related functions.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ src/Internal.hs view
@@ -0,0 +1,5 @@+module Internal ((%^%)) where+import Data.Complex ( Complex (..) )++(%^%) :: Complex Double -> Int -> Complex Double+(%^%) z p = z ^^ p
+ src/Math/Eisenstein.hs view
@@ -0,0 +1,119 @@+module Math.Eisenstein+ ( lambda,+ eisensteinE2,+ eisensteinE4,+ eisensteinE6,+ kleinJ,+ kleinJinv,+ modularDiscriminant,+ agm,+ etaDedekind,+ jtheta1DashDashDash0+ ) where+import Data.Complex ( Complex(..) )+import Internal ( (%^%) )+import Math.EllipticIntegrals ( ellipticF', ellipticE' )+import Math.JacobiTheta ( jtheta2, jtheta3, jtheta4 )+++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 = -2 * etaDedekind tau %^% 3 * eisensteinE2 tau ++
+ src/Math/Weierstrass.hs view
@@ -0,0 +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
+ tests/Approx.hs view
@@ -0,0 +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)
+ tests/Main.hs view
@@ -0,0 +1,144 @@+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++ ]
+ weierstrass-functions.cabal view
@@ -0,0 +1,51 @@+name: weierstrass-functions+version: 0.1.0.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.0.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