packages feed

jacobi-theta (empty) → 0.1.0.0

raw patch · 8 files changed

+312/−0 lines, 8 filesdep +basedep +jacobi-thetadep +tastysetup-changed

Dependencies added: base, jacobi-theta, tasty, tasty-hunit

Files

+ CHANGELOG.md view
@@ -0,0 +1,3 @@+## 0.1.0.0 - 2023-02-17++* Initial 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,3 @@+# jacobi-theta++Evaluation of the Jacobi theta functions.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ jacobi-theta.cabal view
@@ -0,0 +1,45 @@+name:                jacobi-theta+version:             0.1.0.0+synopsis:            Jacobi Theta Functions+description:         Evaluation of the Jacobi theta functions.+homepage:            https://github.com/stla/jacobi-theta#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+cabal-version:       >=1.10+extra-source-files:  README.md+                     CHANGELOG.md++library+  hs-source-dirs:      src+  exposed-modules:     Math.JacobiTheta+  build-depends:       base >= 4.7 && < 5+  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+                      , jacobi-theta+  Default-Language:     Haskell2010++source-repository head+  type:     git+  location: https://github.com/stla/jacobi-theta
+ src/Math/JacobiTheta.hs view
@@ -0,0 +1,147 @@+module Math.JacobiTheta+  (+    jtheta1Dash,+    jtheta1,+    jtheta2,+    jtheta3,+    jtheta4+  )+  where+import Data.Complex++type Cplx = Complex Double++i_ :: Cplx+i_ = 0.0 :+ 1.0++machinePrecision :: Double+machinePrecision = 2**(-52)++areClose :: Cplx -> Cplx -> Bool+areClose z1 z2 = magnitude (z1 - z2) < epsilon * h+  where+    epsilon = 2.0 * machinePrecision+    magn2 = magnitude z2+    h = if magn2 < epsilon then 1.0 else max (magnitude z1) magn2++square :: Cplx -> Cplx+square z = z * z++jtheta1Alt1 :: Cplx -> Cplx -> Cplx+jtheta1Alt1 z q =+  go 0 (0.0 :+ 0.0) 1.0 (1.0 / qsq) 1.0+  where +    qsq = q * q+    go :: Int -> Cplx -> Cplx -> Cplx -> Cplx -> Cplx+    go n out alt q_2n q_n_np1 +      | n > 3000 = error "Reached 3000 iterations."+      | areClose out outnew = 2.0 * sqrt (sqrt q) * out+      | otherwise = go (n + 1) outnew (-alt) q_2np1 q_np1_np2+        where+          q_2np1 = q_2n * qsq+          q_np1_np2 = q_n_np1 * q_2np1+          n' = fromIntegral n +          k = 2.0 * n' + 1.0+          outnew = out + alt * q_np1_np2 * sin (k * z) ++-- jtheta1(z, tau) = jtheta1Alt2 (z/pi) (-i_ * tau/pi)+jtheta1Alt2 :: Cplx -> Cplx -> Cplx+jtheta1Alt2 z' t' = +  let nm = round (0.5 - realPart z') in+  let np = nm + 1 in+  go nm np (0.0 :+ 0.0) (if even np then (-1, 1) else (1, -1)) +  where+    go :: Int -> Int -> Cplx -> (Cplx, Cplx) -> Cplx+    go nminus nplus series (altm, altp)+      | nplus - nminus > 3000 = error "Reached 3000 iterations."+      | (nplus - nminus > 2) && areClose series newseries = +          series / sqrt (pi * t')+      | otherwise = go (nminus - 1) (nplus + 1) newseries (-altm, -altp)+        where +          nminus' = fromIntegral nminus+          nplus' = fromIntegral nplus+          termm = altm * exp (- square (nminus' - 0.5 + z') / t')+          termp = altp * exp (- square (nplus' - 0.5 + z') / t')+          newseries = series + termm + termp++falpha :: Cplx -> Cplx -> Cplx+falpha z tau = +  sqrt (-i_ * tau) * exp (i_ / tau * z * z / pi)++jtheta1Alt :: Cplx -> Cplx -> Cplx+jtheta1Alt z tau = +  if imagPart tau < 1.3 +    then+      let w = pi * tau in +      i_ * jtheta1Alt2 (z / w) (i_ / w) / falpha z tau+    else+      i_ * jtheta1Alt1 (z / tau) (exp (-i_ * pi / tau)) / falpha z tau++tauFromQ :: Cplx -> Cplx+tauFromQ q = -i_ * log q / pi++checkQ :: Cplx -> Cplx+checkQ q+  | magnitude q >= 1 = +    error "The modulus of the nome must be smaller than one."+  | imagPart q == 0 && realPart q <= 0 = +    error "If the nome is real, it must be positive."+  | otherwise = q++getTauFromQ :: Cplx -> Cplx+getTauFromQ = tauFromQ . checkQ++expM :: Cplx -> Cplx -> Cplx+expM z tau = exp (i_ * (z + tau * pi/4))++-- | First Jacobi theta function+jtheta1 :: +     Cplx -- ^ z+  -> Cplx -- ^ q, the nome+  -> Cplx+jtheta1 z q = jtheta1Alt z (getTauFromQ q)++-- | Second Jacobi theta function+jtheta2 :: +     Cplx -- ^ z+  -> Cplx -- ^ q, the nome+  -> Cplx+jtheta2 z = jtheta1 (z + pi/2)++-- | Third Jacobi theta function+jtheta3 :: +     Cplx -- ^ z+  -> Cplx -- ^ q, the nome+  -> Cplx+jtheta3 z q = jtheta2 (z - pi/2 * tau) q * expM (-z) tau+  where+    tau = tauFromQ q++-- | Fourth Jacobi theta function+jtheta4 :: +     Cplx -- ^ z+  -> Cplx -- ^ q, the nome+  -> Cplx+jtheta4 z = jtheta3 (z + pi/2)++-- | Derivative of the first Jacobi theta function+jtheta1Dash :: +     Cplx -- ^ z+  -> Cplx -- ^ q, the nome+  -> Cplx+jtheta1Dash z q = +  go 0 (0.0 :+ 0.0) 1.0 (1.0 / qsq) 1.0+  where +    q' = checkQ q+    qsq = q' * q'+    go :: Int -> Cplx -> Cplx -> Cplx -> Cplx -> Cplx+    go n out alt q_2n q_n_np1 +      | n > 3000 = error "Reached 3000 iterations."+      | areClose out outnew = 2.0 * sqrt (sqrt q) * out+      | otherwise = go (n + 1) outnew (-alt) q_2np1 q_np1_np2+        where+          q_2np1 = q_2n * qsq+          q_np1_np2 = q_n_np1 * q_2np1+          n' = fromIntegral n +          k = 2.0 * n' + 1.0+          outnew = out + k * alt * q_np1_np2 * cos (k * z)
+ tests/Approx.hs view
@@ -0,0 +1,8 @@+module Approx where+import Data.Complex++approx0 :: Int -> Double -> Double+approx0 n x = fromInteger (round $ x * (10^n)) / (10.0^^n)++approx :: Int -> Complex Double -> Complex Double+approx n z = approx0 n (realPart z) :+ approx0 n (imagPart z)
+ tests/Main.hs view
@@ -0,0 +1,74 @@+module Main where+import           Approx+import           Data.Complex+import           Math.JacobiTheta+import           Test.Tasty       (defaultMain, testGroup)+import           Test.Tasty.HUnit (assertEqual, testCase)++i_ :: Complex Double+i_ = 0.0 :+ 1.0++q :: Complex Double +q = exp (-pi)++q' :: Complex Double +q' = exp (-pi/100)++q'' :: Complex Double +q'' = exp (i_ * pi * tau)+  where+    tau = 2.0 :+ 2.0++main :: IO ()+main = defaultMain $+  testGroup "Tests"+  [ testCase "a jtheta1 value" $ do+      let expected = 1.1816128551455719 :+ 0.59589712760417439+          obtained = jtheta1 (1 :+ 1) q+      assertEqual ""+        (approx 10 obtained)+        (approx 10 expected),++    testCase "another jtheta1 value" $ do+      let expected = 0.0284051242069853 :+ 0.0+          obtained = jtheta1 2 q'+      assertEqual ""+        (approx 10 obtained)+        (approx 10 expected),++    testCase "yet another jtheta1 value" $ do+      let expected = 0.539843563932874 :+ 0.26400643871132+          obtained = jtheta1 (1 :+ 1) q''+      assertEqual ""+        (approx 10 obtained)+        (approx 10 expected),++    testCase "a jtheta2 value" $ do+      let expected = 0.74328632006610539 :+ (-0.904159309718008)+          obtained = jtheta2 (1 :+ 1) q+      assertEqual ""+        (approx 10 obtained)+        (approx 10 expected),++    testCase "a jtheta3 value" $ do+      let expected = 0.86456184935441778 :+ (-0.28488586703507289)+          obtained = jtheta3 (1 :+ 1) q+      assertEqual ""+        (approx 10 obtained)+        (approx 10 expected),++    testCase "a jtheta4 value" $ do+      let expected = 1.1351891564632007 :+ 0.28517396444192509+          obtained = jtheta4 (1 :+ 1) q+      assertEqual ""+        (approx 10 obtained)+        (approx 10 expected),++    testCase "a jtheta1Dash value" $ do+      let expected = 0.81117649363854416 :+ (-0.89452803853474627)+          obtained = jtheta1Dash (1 :+ 1) q+      assertEqual ""+        (approx 10 obtained)+        (approx 10 expected)++  ]