Open stla opened 2 years ago
Since q = exp(iπτ) = 0.999999685... in this case and the theta functions are only defined for |q| < 1, it isn't surprising that there will be some convergence problems if one approaches this boundary too closely. Of course, Wolfram can do it, but maybe they need a lot of iterations there too.
@simonp0420 I've ported the new Fortran implementation by @fremling to Haskell and with this implementation the result is immediately obtained! That's impressive. I think I will port it to R/C++ now, and if all unit tests succeed, maybe in Julia later.
Could you provide a link to the Fortran implementation so I could take a look?
Here. He uses the modular transformations to bring $\tau$ to a certain region where the evaluation works.
@simonp0420 I ported the code to R. It works, and better than the previous one!
I added some comments to the gist.
I found a value for which the 3000 iterations are attained when calculating
jtheta2
: $$z=0 \quad \text{and} \quad \tau = 0.7792256 +10^{-7}i.$$This is really a edge case: with
tau = 0.78 + 1e-7
, this works.There's a formula linking
theta(z, tau)
andtheta(z/tau, -1/tau)
, and one can get the result with this formula. One can see the presence of the edge case once again, because the imaginary part of-1/tau
is1.646924e-07
.@simonp0420 do you have any comment ?
See this issue for a discussion about this problem, using a different implementation (the one I used before @simonp0420 proposed the new implementation).