module Math.JacobiTheta
(
jtheta1,
jtheta2,
jtheta3,
jtheta4,
jthetaAB,
jtheta1Dash0,
jtheta1Dash
)
where
import Data.Complex ( imagPart, magnitude, realPart, Complex(..) )
type Cplx = Complex Double
i_ :: Cplx
i_ :: Complex Double
i_ = Double
0.0 forall a. a -> a -> Complex a
:+ Double
1.0
machinePrecision :: Double
machinePrecision :: Double
machinePrecision = Double
2forall a. Floating a => a -> a -> a
**(-Double
52)
areClose :: Cplx -> Cplx -> Bool
areClose :: Complex Double -> Complex Double -> Bool
areClose Complex Double
z1 Complex Double
z2 = forall a. RealFloat a => Complex a -> a
magnitude (Complex Double
z1 forall a. Num a => a -> a -> a
- Complex Double
z2) forall a. Ord a => a -> a -> Bool
< Double
epsilon forall a. Num a => a -> a -> a
* Double
h
where
epsilon :: Double
epsilon = Double
2.0 forall a. Num a => a -> a -> a
* Double
machinePrecision
magn2 :: Double
magn2 = forall a. RealFloat a => Complex a -> a
magnitude Complex Double
z2
h :: Double
h = if Double
magn2 forall a. Ord a => a -> a -> Bool
< Double
epsilon then Double
1.0 else forall a. Ord a => a -> a -> a
max (forall a. RealFloat a => Complex a -> a
magnitude Complex Double
z1) Double
magn2
modulo :: Double -> Int -> Double
modulo :: Double -> Int -> Double
modulo Double
a Int
p =
let p' :: Double
p' = forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
p
in
if Double
a forall a. Ord a => a -> a -> Bool
> Double
0
then Double
a forall a. Num a => a -> a -> a
- forall a b. (Integral a, Num b) => a -> b
fromIntegral(Int
p forall a. Num a => a -> a -> a
* forall a b. (RealFrac a, Integral b) => a -> b
floor(Double
aforall a. Fractional a => a -> a -> a
/Double
p'))
else Double
a forall a. Num a => a -> a -> a
- forall a b. (Integral a, Num b) => a -> b
fromIntegral(Int
p forall a. Num a => a -> a -> a
* forall a b. (RealFrac a, Integral b) => a -> b
ceiling(Double
aforall a. Fractional a => a -> a -> a
/Double
p'))
dologtheta4 :: Cplx -> Cplx -> Int -> Int -> Cplx
dologtheta4 :: Complex Double -> Complex Double -> Int -> Int -> Complex Double
dologtheta4 Complex Double
z Complex Double
tau Int
passes Int
maxiter =
Complex Double -> Complex Double -> Int -> Int -> Complex Double
dologtheta3 (Complex Double
z forall a. Num a => a -> a -> a
+ Complex Double
0.5) Complex Double
tau (Int
passesforall a. Num a => a -> a -> a
+Int
1) Int
maxiter
dologtheta3 :: Cplx -> Cplx -> Int -> Int -> Cplx
dologtheta3 :: Complex Double -> Complex Double -> Int -> Int -> Complex Double
dologtheta3 Complex Double
z Complex Double
tau Int
passes Int
maxiterloc
| forall a. Complex a -> a
realPart Complex Double
tau2 forall a. Ord a => a -> a -> Bool
> Double
0.6 = Complex Double -> Complex Double -> Int -> Int -> Complex Double
dologtheta4 Complex Double
z (Complex Double
tau2 forall a. Num a => a -> a -> a
- Complex Double
1) (Int
passes forall a. Num a => a -> a -> a
+ Int
1) Int
maxiterloc
| forall a. Complex a -> a
realPart Complex Double
tau2 forall a. Ord a => a -> a -> Bool
< -Double
0.6 = Complex Double -> Complex Double -> Int -> Int -> Complex Double
dologtheta4 Complex Double
z (Complex Double
tau2 forall a. Num a => a -> a -> a
+ Complex Double
1) (Int
passes forall a. Num a => a -> a -> a
+ Int
1) Int
maxiterloc
| forall a. RealFloat a => Complex a -> a
magnitude Complex Double
tau2 forall a. Ord a => a -> a -> Bool
< Double
0.98 Bool -> Bool -> Bool
&& forall a. Complex a -> a
imagPart Complex Double
tau2 forall a. Ord a => a -> a -> Bool
< Double
0.98 =
Complex Double
i_ forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* Complex Double
tauprime forall a. Num a => a -> a -> a
* Complex Double
z forall a. Num a => a -> a -> a
* Complex Double
z
forall a. Num a => a -> a -> a
+ Complex Double -> Complex Double -> Int -> Int -> Complex Double
dologtheta3 (Complex Double
z forall a. Num a => a -> a -> a
* Complex Double
tauprime) Complex Double
tauprime (Int
passes forall a. Num a => a -> a -> a
+ Int
1) Int
maxiterloc
forall a. Num a => a -> a -> a
- forall a. Floating a => a -> a
log(forall a. Floating a => a -> a
sqrt Complex Double
tau2 forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
sqrt Complex Double
i_)
| Bool
otherwise = Complex Double -> Complex Double -> Int -> Int -> Complex Double
argtheta3 Complex Double
z Complex Double
tau2 Int
0 Int
maxiterloc
where
rPtau :: Double
rPtau = forall a. Complex a -> a
realPart Complex Double
tau
rPtau2 :: Double
rPtau2 = if Double
rPtau forall a. Ord a => a -> a -> Bool
> Double
0
then Double -> Int -> Double
modulo (Double
rPtau forall a. Num a => a -> a -> a
+ Double
1) Int
2 forall a. Num a => a -> a -> a
- Double
1
else Double -> Int -> Double
modulo (Double
rPtau forall a. Num a => a -> a -> a
- Double
1) Int
2 forall a. Num a => a -> a -> a
+ Double
1
tau2 :: Complex Double
tau2 = Double
rPtau2 forall a. a -> a -> Complex a
:+ forall a. Complex a -> a
imagPart Complex Double
tau
tauprime :: Complex Double
tauprime = -Complex Double
1 forall a. Fractional a => a -> a -> a
/ Complex Double
tau2
argtheta3 :: Cplx -> Cplx -> Int -> Int -> Cplx
argtheta3 :: Complex Double -> Complex Double -> Int -> Int -> Complex Double
argtheta3 Complex Double
z Complex Double
tau Int
passes Int
maxiterloc
| Int
passes forall a. Ord a => a -> a -> Bool
> Int
maxiterloc = forall a. HasCallStack => [Char] -> a
error [Char]
"Reached maximal iteration."
| Double
iPz forall a. Ord a => a -> a -> Bool
< -Double
iPtau forall a. Fractional a => a -> a -> a
/ Double
2 = Complex Double -> Complex Double -> Int -> Int -> Complex Double
argtheta3 (-Complex Double
zuse) Complex Double
tau (Int
passes forall a. Num a => a -> a -> a
+ Int
1) Int
maxiterloc
| Double
iPz forall a. Ord a => a -> a -> Bool
>= Double
iPtau forall a. Fractional a => a -> a -> a
/ Double
2 =
-Complex Double
2 forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* Complex Double
quotient forall a. Num a => a -> a -> a
* Complex Double
i_ forall a. Num a => a -> a -> a
* Complex Double
zmin
forall a. Num a => a -> a -> a
+ Complex Double -> Complex Double -> Int -> Int -> Complex Double
argtheta3 Complex Double
zmin Complex Double
tau (Int
passes forall a. Num a => a -> a -> a
+ Int
1) Int
maxiterloc
forall a. Num a => a -> a -> a
- Complex Double
i_ forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* Complex Double
tau forall a. Num a => a -> a -> a
* Complex Double
quotient forall a. Num a => a -> a -> a
* Complex Double
quotient
| Bool
otherwise = Complex Double -> Complex Double -> Complex Double
calctheta3 Complex Double
zuse Complex Double
tau
where
iPz :: Double
iPz = forall a. Complex a -> a
imagPart Complex Double
z
iPtau :: Double
iPtau = forall a. Complex a -> a
imagPart Complex Double
tau
zuse :: Complex Double
zuse = Double -> Int -> Double
modulo (forall a. Complex a -> a
realPart Complex Double
z) Int
1 forall a. a -> a -> Complex a
:+ Double
iPz
quotient :: Complex Double
quotient = Int -> Complex Double
fromInt forall a b. (a -> b) -> a -> b
$ forall a b. (RealFrac a, Integral b) => a -> b
floor(Double
iPz forall a. Fractional a => a -> a -> a
/ Double
iPtau forall a. Num a => a -> a -> a
+ Double
0.5)
zmin :: Complex Double
zmin = Complex Double
zuse forall a. Num a => a -> a -> a
- Complex Double
tau forall a. Num a => a -> a -> a
* Complex Double
quotient
fromInt :: Int -> Cplx
fromInt :: Int -> Complex Double
fromInt = forall a b. (Integral a, Num b) => a -> b
fromIntegral
calctheta3 :: Cplx -> Cplx -> Cplx
calctheta3 :: Complex Double -> Complex Double -> Complex Double
calctheta3 Complex Double
z Complex Double
tau =
Int -> Complex Double -> Complex Double
go Int
1 Complex Double
1
where
qw :: Int -> Cplx
qw :: Int -> Complex Double
qw Int
n = forall a. Floating a => a -> a
exp(Complex Double
inpi forall a. Num a => a -> a -> a
* (Complex Double
taun forall a. Num a => a -> a -> a
+ Complex Double
2 forall a. Num a => a -> a -> a
* Complex Double
z)) forall a. Num a => a -> a -> a
+ forall a. Floating a => a -> a
exp(Complex Double
inpi forall a. Num a => a -> a -> a
* (Complex Double
taun forall a. Num a => a -> a -> a
- Complex Double
2 forall a. Num a => a -> a -> a
* Complex Double
z))
where
n' :: Complex Double
n' = forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
inpi :: Complex Double
inpi = Complex Double
i_ forall a. Num a => a -> a -> a
* Complex Double
n' forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi
taun :: Complex Double
taun = Complex Double
n' forall a. Num a => a -> a -> a
* Complex Double
tau
go :: Int -> Complex Double -> Complex Double
go Int
n Complex Double
res
| forall a. RealFloat a => a -> Bool
isNaN Double
modulus = forall a. HasCallStack => [Char] -> a
error [Char]
"NaN has occured in the summation."
| forall a. RealFloat a => a -> Bool
isInfinite Double
modulus = forall a. HasCallStack => [Char] -> a
error [Char]
"Infinity reached in the summation."
| Int
n forall a. Ord a => a -> a -> Bool
>= Int
3 Bool -> Bool -> Bool
&& Complex Double -> Complex Double -> Bool
areClose Complex Double
res Complex Double
resnew = forall a. Floating a => a -> a
log Complex Double
res
| Bool
otherwise = Int -> Complex Double -> Complex Double
go (Int
n forall a. Num a => a -> a -> a
+ Int
1) Complex Double
resnew
where
modulus :: Double
modulus = forall a. RealFloat a => Complex a -> a
magnitude Complex Double
res
resnew :: Complex Double
resnew = Complex Double
res forall a. Num a => a -> a -> a
+ Int -> Complex Double
qw Int
n
tauFromQ :: Cplx -> Cplx
tauFromQ :: Complex Double -> Complex Double
tauFromQ Complex Double
q = -Complex Double
i_ forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
log Complex Double
q forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a
pi
checkQ :: Cplx -> Cplx
checkQ :: Complex Double -> Complex Double
checkQ Complex Double
q
| forall a. RealFloat a => Complex a -> a
magnitude Complex Double
q forall a. Ord a => a -> a -> Bool
>= Double
1 =
forall a. HasCallStack => [Char] -> a
error [Char]
"The modulus of the nome must be smaller than one."
| forall a. Complex a -> a
imagPart Complex Double
q forall a. Eq a => a -> a -> Bool
== Double
0 Bool -> Bool -> Bool
&& forall a. Complex a -> a
realPart Complex Double
q forall a. Ord a => a -> a -> Bool
<= Double
0 =
forall a. HasCallStack => [Char] -> a
error [Char]
"If the nome is real, it must be positive."
| Bool
otherwise = Complex Double
q
getTauFromQ :: Cplx -> Cplx
getTauFromQ :: Complex Double -> Complex Double
getTauFromQ = Complex Double -> Complex Double
tauFromQ forall b c a. (b -> c) -> (a -> b) -> a -> c
. Complex Double -> Complex Double
checkQ
funM :: Cplx -> Cplx -> Cplx
funM :: Complex Double -> Complex Double -> Complex Double
funM Complex Double
z Complex Double
tau = Complex Double
i_ forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* (Complex Double
z forall a. Num a => a -> a -> a
+ Complex Double
tauforall a. Fractional a => a -> a -> a
/Complex Double
4)
ljtheta1 :: Cplx -> Cplx -> Cplx
ljtheta1 :: Complex Double -> Complex Double -> Complex Double
ljtheta1 Complex Double
z Complex Double
tau = Complex Double -> Complex Double -> Complex Double
ljtheta2 (Complex Double
z forall a. Num a => a -> a -> a
- Complex Double
0.5) Complex Double
tau
jtheta1 ::
Complex Double
-> Complex Double
-> Complex Double
jtheta1 :: Complex Double -> Complex Double -> Complex Double
jtheta1 Complex Double
z Complex Double
q = forall a. Floating a => a -> a
exp(Complex Double -> Complex Double -> Complex Double
ljtheta1 (Complex Double
zforall a. Fractional a => a -> a -> a
/forall a. Floating a => a
pi) Complex Double
tau)
where
tau :: Complex Double
tau = Complex Double -> Complex Double
getTauFromQ Complex Double
q
ljtheta2 :: Cplx -> Cplx -> Cplx
ljtheta2 :: Complex Double -> Complex Double -> Complex Double
ljtheta2 Complex Double
z Complex Double
tau =
Complex Double -> Complex Double -> Complex Double
funM Complex Double
z Complex Double
tau forall a. Num a => a -> a -> a
+ Complex Double -> Complex Double -> Int -> Int -> Complex Double
dologtheta3 (Complex Double
z forall a. Num a => a -> a -> a
+ Complex Double
0.5 forall a. Num a => a -> a -> a
* Complex Double
tau) Complex Double
tau Int
0 Int
1000
jtheta2 ::
Complex Double
-> Complex Double
-> Complex Double
jtheta2 :: Complex Double -> Complex Double -> Complex Double
jtheta2 Complex Double
z Complex Double
q = forall a. Floating a => a -> a
exp(Complex Double -> Complex Double -> Complex Double
ljtheta2 (Complex Double
zforall a. Fractional a => a -> a -> a
/forall a. Floating a => a
pi) Complex Double
tau)
where
tau :: Complex Double
tau = Complex Double -> Complex Double
getTauFromQ Complex Double
q
jtheta3 ::
Complex Double
-> Complex Double
-> Complex Double
jtheta3 :: Complex Double -> Complex Double -> Complex Double
jtheta3 Complex Double
z Complex Double
q = forall a. Floating a => a -> a
exp(Complex Double -> Complex Double -> Int -> Int -> Complex Double
dologtheta3 (Complex Double
zforall a. Fractional a => a -> a -> a
/forall a. Floating a => a
pi) Complex Double
tau Int
0 Int
1000)
where
tau :: Complex Double
tau = Complex Double -> Complex Double
getTauFromQ Complex Double
q
jtheta4 ::
Complex Double
-> Complex Double
-> Complex Double
jtheta4 :: Complex Double -> Complex Double -> Complex Double
jtheta4 Complex Double
z Complex Double
q = forall a. Floating a => a -> a
exp(Complex Double -> Complex Double -> Int -> Int -> Complex Double
dologtheta4 (Complex Double
zforall a. Fractional a => a -> a -> a
/forall a. Floating a => a
pi) Complex Double
tau Int
0 Int
1000)
where
tau :: Complex Double
tau = Complex Double -> Complex Double
getTauFromQ Complex Double
q
jthetaAB' ::
Complex Double
-> Complex Double
-> Complex Double
-> Complex Double
-> Complex Double
jthetaAB' :: Complex Double
-> Complex Double
-> Complex Double
-> Complex Double
-> Complex Double
jthetaAB' Complex Double
a Complex Double
b Complex Double
z Complex Double
tau = Complex Double
c forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
exp(Complex Double -> Complex Double -> Int -> Int -> Complex Double
dologtheta3 (Complex Double
alphaforall a. Num a => a -> a -> a
+Complex Double
beta) Complex Double
tau Int
0 Int
1000)
where
alpha :: Complex Double
alpha = Complex Double
a forall a. Num a => a -> a -> a
* Complex Double
tau
beta :: Complex Double
beta = Complex Double
zforall a. Fractional a => a -> a -> a
/forall a. Floating a => a
pi forall a. Num a => a -> a -> a
+ Complex Double
b
c :: Complex Double
c = forall a. Floating a => a -> a
exp(Complex Double
i_ forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* Complex Double
a forall a. Num a => a -> a -> a
* (Complex Double
alpha forall a. Num a => a -> a -> a
+ Complex Double
2forall a. Num a => a -> a -> a
*Complex Double
beta))
jthetaAB ::
Complex Double
-> Complex Double
-> Complex Double
-> Complex Double
-> Complex Double
jthetaAB :: Complex Double
-> Complex Double
-> Complex Double
-> Complex Double
-> Complex Double
jthetaAB Complex Double
a Complex Double
b Complex Double
z Complex Double
q = Complex Double
c forall a. Num a => a -> a -> a
* Complex Double -> Complex Double -> Complex Double
jtheta3 (Complex Double
alpha forall a. Num a => a -> a -> a
+ Complex Double
beta) Complex Double
q
where
tau :: Complex Double
tau = Complex Double -> Complex Double
getTauFromQ Complex Double
q
alpha :: Complex Double
alpha = forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* Complex Double
a forall a. Num a => a -> a -> a
* Complex Double
tau
beta :: Complex Double
beta = Complex Double
z forall a. Num a => a -> a -> a
+ forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* Complex Double
b
c :: Complex Double
c = forall a. Floating a => a -> a
exp(Complex Double
i_ forall a. Num a => a -> a -> a
* Complex Double
a forall a. Num a => a -> a -> a
* (Complex Double
alpha forall a. Num a => a -> a -> a
+ Complex Double
2forall a. Num a => a -> a -> a
*Complex Double
beta))
jtheta1Dash0 ::
Complex Double
-> Complex Double
jtheta1Dash0 :: Complex Double -> Complex Double
jtheta1Dash0 Complex Double
q =
-Complex Double
2 forall a. Num a => a -> a -> a
* Complex Double
i_ forall a. Num a => a -> a -> a
* Complex Double
jab forall a. Num a => a -> a -> a
* Complex Double
jab forall a. Num a => a -> a -> a
* Complex Double
jab
where
tau :: Complex Double
tau = Complex Double -> Complex Double
getTauFromQ Complex Double
q
jab :: Complex Double
jab = Complex Double
-> Complex Double
-> Complex Double
-> Complex Double
-> Complex Double
jthetaAB' (Complex Double
1forall a. Fractional a => a -> a -> a
/Complex Double
6) Complex Double
0.5 Complex Double
0 (Complex Double
3forall a. Num a => a -> a -> a
*Complex Double
tau)
jtheta1Dash ::
Complex Double
-> Complex Double
-> Complex Double
jtheta1Dash :: Complex Double -> Complex Double -> Complex Double
jtheta1Dash Complex Double
z Complex Double
q =
Int
-> Complex Double
-> Complex Double
-> Complex Double
-> Complex Double
-> Complex Double
go Int
0 (Double
0.0 forall a. a -> a -> Complex a
:+ Double
0.0) Complex Double
1.0 (Complex Double
1.0 forall a. Fractional a => a -> a -> a
/ Complex Double
qsq) Complex Double
1.0
where
q' :: Complex Double
q' = Complex Double -> Complex Double
checkQ Complex Double
q
qsq :: Complex Double
qsq = Complex Double
q' forall a. Num a => a -> a -> a
* Complex Double
q'
go :: Int -> Cplx -> Cplx -> Cplx -> Cplx -> Cplx
go :: Int
-> Complex Double
-> Complex Double
-> Complex Double
-> Complex Double
-> Complex Double
go Int
n Complex Double
out Complex Double
alt Complex Double
q_2n Complex Double
q_n_np1
| Int
n forall a. Ord a => a -> a -> Bool
> Int
3000 = forall a. HasCallStack => [Char] -> a
error [Char]
"Reached 3000 iterations."
| Complex Double -> Complex Double -> Bool
areClose Complex Double
out Complex Double
outnew = Complex Double
2.0 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sqrt (forall a. Floating a => a -> a
sqrt Complex Double
q) forall a. Num a => a -> a -> a
* Complex Double
out
| Bool
otherwise = Int
-> Complex Double
-> Complex Double
-> Complex Double
-> Complex Double
-> Complex Double
go (Int
n forall a. Num a => a -> a -> a
+ Int
1) Complex Double
outnew (-Complex Double
alt) Complex Double
q_2np1 Complex Double
q_np1_np2
where
q_2np1 :: Complex Double
q_2np1 = Complex Double
q_2n forall a. Num a => a -> a -> a
* Complex Double
qsq
q_np1_np2 :: Complex Double
q_np1_np2 = Complex Double
q_n_np1 forall a. Num a => a -> a -> a
* Complex Double
q_2np1
n' :: Complex Double
n' = forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
k :: Complex Double
k = Complex Double
2.0 forall a. Num a => a -> a -> a
* Complex Double
n' forall a. Num a => a -> a -> a
+ Complex Double
1.0
outnew :: Complex Double
outnew = Complex Double
out forall a. Num a => a -> a -> a
+ Complex Double
k forall a. Num a => a -> a -> a
* Complex Double
alt forall a. Num a => a -> a -> a
* Complex Double
q_np1_np2 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos (Complex Double
k forall a. Num a => a -> a -> a
* Complex Double
z)