diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,3 @@ +## 0.1.0.0 - 2023-02-17 + +* Initial release. diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,3 @@ +# jacobi-theta + +Evaluation of the Jacobi theta functions. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/jacobi-theta.cabal b/jacobi-theta.cabal new file mode 100644 --- /dev/null +++ b/jacobi-theta.cabal @@ -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 diff --git a/src/Math/JacobiTheta.hs b/src/Math/JacobiTheta.hs new file mode 100644 --- /dev/null +++ b/src/Math/JacobiTheta.hs @@ -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) diff --git a/tests/Approx.hs b/tests/Approx.hs new file mode 100644 --- /dev/null +++ b/tests/Approx.hs @@ -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) diff --git a/tests/Main.hs b/tests/Main.hs new file mode 100644 --- /dev/null +++ b/tests/Main.hs @@ -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) + + ]