jacobi-theta 0.1.2.0 → 0.2.0.0
raw patch · 8 files changed
+384/−333 lines, 8 filessetup-changedPVP ok
version bump matches the API change (PVP)
API changes (from Hackage documentation)
Files
- CHANGELOG.md +23/−18
- LICENSE +30/−30
- README.md +7/−7
- Setup.hs +2/−2
- jacobi-theta.cabal +44/−44
- src/Math/JacobiTheta.hs +175/−150
- tests/Approx.hs +8/−8
- tests/Main.hs +95/−74
CHANGELOG.md view
@@ -1,18 +1,23 @@-## 0.1.2.0 - 2023-03-02--Some values of the Jacobi theta functions were wrong.---## 0.1.1.1 - 2023-03-01--Removed a useless type alias.---## 0.1.1.0 - 2023-02-25--Some values of the Jacobi theta functions were wrong.---## 0.1.0.0 - 2023-02-17--Initial release.+## 0.2.0.0 - 2023-10-16 + +New, better implementation of the Jacobi theta functions. + + +## 0.1.2.0 - 2023-03-02 + +Some values of the Jacobi theta functions were wrong. + + +## 0.1.1.1 - 2023-03-01 + +Removed a useless type alias. + + +## 0.1.1.0 - 2023-02-25 + +Some values of the Jacobi theta functions were wrong. + + +## 0.1.0.0 - 2023-02-17 + +Initial 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,7 +1,7 @@-# jacobi-theta--<!-- badges: start -->-[](https://github.com/stla/jacobi-theta/actions/workflows/GithubAction.yml)-<!-- badges: end -->--Evaluation of the Jacobi theta functions.+# jacobi-theta + +<!-- badges: start --> +[](https://github.com/stla/jacobi-theta/actions/workflows/GithubAction.yml) +<!-- badges: end --> + +Evaluation of the Jacobi theta functions.
Setup.hs view
@@ -1,2 +1,2 @@-import Distribution.Simple-main = defaultMain+import Distribution.Simple +main = defaultMain
jacobi-theta.cabal view
@@ -1,45 +1,45 @@-name: jacobi-theta-version: 0.1.2.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+name: jacobi-theta +version: 0.2.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
@@ -1,150 +1,175 @@-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- jtheta1Alt1 z (exp (i_ * pi * tau))--- let w = pi * tau in --- i_ * jtheta1Alt2 (z / w) (i_ / w) / falpha z tau- else- let t = - i_ * tau / pi in- jtheta1Alt2 (z / pi) t--- 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 :: - Complex Double -- ^ z- -> Complex Double -- ^ q, the nome- -> Complex Double-jtheta1 z q = jtheta1Alt z (getTauFromQ q)---- | Second Jacobi theta function-jtheta2 :: - Complex Double -- ^ z- -> Complex Double -- ^ q, the nome- -> Complex Double-jtheta2 z = jtheta1 (z + pi/2)---- | Third Jacobi theta function-jtheta3 :: - Complex Double -- ^ z- -> Complex Double -- ^ q, the nome- -> Complex Double-jtheta3 z q = jtheta2 (z - pi/2 * tau) q * expM (-z) tau- where- tau = tauFromQ q---- | Fourth Jacobi theta function-jtheta4 :: - Complex Double -- ^ z- -> Complex Double -- ^ q, the nome- -> Complex Double-jtheta4 z = jtheta3 (z + pi/2)---- | Derivative of the first Jacobi theta function-jtheta1Dash :: - Complex Double -- ^ z- -> Complex Double -- ^ q, the nome- -> Complex Double-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) +module Math.JacobiTheta + ( + jtheta1, + jtheta2, + jtheta3, + jtheta4, + jtheta1Dash + ) + where +import Data.Complex ( imagPart, magnitude, realPart, 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 + +modulo :: Double -> Int -> Double +modulo a p = + let p' = fromIntegral p + in + if a > 0 + then a - fromIntegral(p * floor(a/p')) + else a - fromIntegral(p * ceiling(a/p')) + +dologtheta4 :: Cplx -> Cplx -> Int -> Int -> Cplx +dologtheta4 z tau passes maxiter = + dologtheta3 (z + 0.5) tau (passes+1) maxiter + +dologtheta3 :: Cplx -> Cplx -> Int -> Int -> Cplx +dologtheta3 z tau passes maxiterloc + | realPart tau2 > 0.6 = dologtheta4 z (tau2 - 1) (passes + 1) maxiterloc + | realPart tau2 < -0.6 = dologtheta4 z (tau2 + 1) (passes + 1) maxiterloc + | magnitude tau2 < 0.98 && imagPart tau2 < 0.98 = + i_ * pi * tauprime * z * z + + dologtheta3 (z * tauprime) tauprime (passes + 1) maxiterloc + - log(sqrt tau2 / sqrt i_) + | otherwise = argtheta3 z tau2 0 maxiterloc + where + rPtau = realPart tau + rPtau2 = if rPtau > 0 + then modulo (rPtau + 1) 2 - 1 + else modulo (rPtau - 1) 2 + 1 + tau2 = rPtau2 :+ imagPart tau + tauprime = -1 / tau2 + +argtheta3 :: Cplx -> Cplx -> Int -> Int -> Cplx +argtheta3 z tau passes maxiterloc + | passes > maxiterloc = error "Reached maximal iteration." + | iPz < -iPtau / 2 = argtheta3 (-zuse) tau (passes + 1) maxiterloc + | iPz >= iPtau / 2 = + -2 * pi * quotient * i_ * zmin + + argtheta3 zmin tau (passes + 1) maxiterloc + - i_ * pi * tau * quotient * quotient + | otherwise = calctheta3 zuse tau + where + iPz = imagPart z + iPtau = imagPart tau + zuse = modulo (realPart z) 1 :+ iPz + quotient = fromInt $ floor(iPz / iPtau + 0.5) + zmin = zuse - tau * quotient + fromInt :: Int -> Cplx + fromInt = fromIntegral + +calctheta3 :: Cplx -> Cplx -> Cplx +calctheta3 z tau = + go 1 1 + where + qw :: Int -> Cplx + qw n = exp(inpi * (taun + 2 * z)) + exp(inpi * (taun - 2 * z)) + where + n' = fromIntegral n + inpi = i_ * n' * pi + taun = n' * tau + go n res + | isNaN modulus = error "NaN has occured in the summation." + | isInfinite modulus = error "Infinity reached in the summation." +-- | modulus == 0 = error "Zero has occured in the summation." + | n >= 3 && areClose res resnew = log res + | otherwise = go (n + 1) resnew + where + modulus = magnitude res + resnew = res + qw n + +------------------------------------------------------------------------------- +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 + +funM :: Cplx -> Cplx -> Cplx +funM z tau = i_ * pi * (z + tau/4) + +ljtheta1 :: Cplx -> Cplx -> Cplx +ljtheta1 z tau = ljtheta2 (z - 0.5) tau + +-- | First Jacobi theta function +jtheta1 :: + Complex Double -- ^ z + -> Complex Double -- ^ q, the nome + -> Complex Double +jtheta1 z q = exp(ljtheta1 (z/pi) tau) + where + tau = getTauFromQ q + +ljtheta2 :: Cplx -> Cplx -> Cplx +ljtheta2 z tau = + funM z tau + dologtheta3 (z + 0.5 * tau) tau 0 1000 + +-- | Second Jacobi theta function +jtheta2 :: + Complex Double -- ^ z + -> Complex Double -- ^ q, the nome + -> Complex Double +jtheta2 z q = exp(ljtheta2 (z/pi) tau) + where + tau = getTauFromQ q + +-- | Third Jacobi theta function +jtheta3 :: + Complex Double -- ^ z + -> Complex Double -- ^ q, the nome + -> Complex Double +jtheta3 z q = exp(dologtheta3 (z/pi) tau 0 1000) + where + tau = getTauFromQ q + +-- | Fourth Jacobi theta function +jtheta4 :: + Complex Double -- ^ z + -> Complex Double -- ^ q, the nome + -> Complex Double +jtheta4 z q = exp(dologtheta4 (z/pi) tau 0 1000) + where + tau = getTauFromQ q + +-- | Derivative of the first Jacobi theta function +jtheta1Dash :: + Complex Double -- ^ z + -> Complex Double -- ^ q, the nome + -> Complex Double +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
@@ -1,8 +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)+module Approx where +import Data.Complex ( imagPart, realPart, 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
@@ -1,74 +1,95 @@-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/10)--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 = 1.75929905417707 :+ 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)-- ]+module Main where +import Approx ( approx ) +import Data.Complex ( Complex(..) ) +import Math.JacobiTheta + ( jtheta1, jtheta2, jtheta3, jtheta4, jtheta1Dash ) +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/10) + +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 = 1.75929905417707 :+ 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), + + testCase "Jacobi identity" $ do + let q''' = 0.556 :+ 0.283 + theta2 = jtheta2 0 q''' + theta3 = jtheta3 0 q''' + theta4 = jtheta4 0 q''' + expected = theta3**4 + obtained = theta2**4 + theta4**4 + assertEqual "" + (approx 10 obtained) + (approx 10 expected), + + testCase "Edge case" $ do + let tau = 0.7792256 :+ 1.0e-7 + q'''' = exp(i_ * pi * tau) + obtained = jtheta2 0 q'''' + expected = 27.746815969548447 :+ 31.241216782108797 + assertEqual "" + (approx 10 obtained) + (approx 10 expected) + + ]