diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -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. diff --git a/LICENSE b/LICENSE --- a/LICENSE +++ b/LICENSE @@ -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. diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -1,7 +1,7 @@ -# jacobi-theta - - -[![Check](https://github.com/stla/jacobi-theta/actions/workflows/GithubAction.yml/badge.svg)](https://github.com/stla/jacobi-theta/actions/workflows/GithubAction.yml) - - -Evaluation of the Jacobi theta functions. +# jacobi-theta + + +[![Check](https://github.com/stla/jacobi-theta/actions/workflows/GithubAction.yml/badge.svg)](https://github.com/stla/jacobi-theta/actions/workflows/GithubAction.yml) + + +Evaluation of the Jacobi theta functions. diff --git a/Setup.hs b/Setup.hs --- a/Setup.hs +++ b/Setup.hs @@ -1,2 +1,2 @@ -import Distribution.Simple -main = defaultMain +import Distribution.Simple +main = defaultMain diff --git a/jacobi-theta.cabal b/jacobi-theta.cabal --- a/jacobi-theta.cabal +++ b/jacobi-theta.cabal @@ -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 diff --git a/src/Math/JacobiTheta.hs b/src/Math/JacobiTheta.hs --- a/src/Math/JacobiTheta.hs +++ b/src/Math/JacobiTheta.hs @@ -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) diff --git a/tests/Approx.hs b/tests/Approx.hs --- a/tests/Approx.hs +++ b/tests/Approx.hs @@ -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) diff --git a/tests/Main.hs b/tests/Main.hs --- a/tests/Main.hs +++ b/tests/Main.hs @@ -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) + + ]