# SiegelTheta provides a wrong value?

Bug introduced in 6.0 or earlier and persists through 11.0.1 or later

I find that SiegelTheta provides a wrong value for high precision. For example, I have the following matrix:

P = {{2 I, I, -3, 3}, {I, 2 I, 0, -3}, {-3, 0, 2 I, I}, {3, -3, I, 2 I}}


and I want to evaluate SiegelTheta[P, {0, 0, 0, 0}]. According to Mathematica, I get values which differ a lot from the correct value as I increase the precision:

N[SiegelTheta[P, {0, 0, 0, 0}], 50]
1.02253493627747251337360506390118996977716349302611 + 0.*10^-51 I

N[SiegelTheta[P, {0, 0, 0, 0}], 60]
1.022534936277472513373605063901189969777163493026107090953727 + 0.*10^-61 I

N[SiegelTheta[P, {0, 0, 0, 0}], 70]
1072.0058459858069455194722642419955529253887586996976403456379723195717 + 0.*10^-68 I

N[SiegelTheta[P, {0, 0, 0, 0}], 100]
4.895826585168781732011692984399098064056816247435531538659849186987471892518027397637334165655318210*10^46 + 0.*10^-54 I


N[SiegelTheta[P, {0, 0, 0, 0}], 50] is $O(1)$ and N[SiegelTheta[P, {0, 0, 0, 0}], 100] is $O(10^{46})$, and the one has very very different order from the other. I can't understand why these two values are very different. (This phenomenon arises only when we use $n\times n$ matrices where $n>3$.)

I want to correctly evaluate SiegelTheta for high precision. I would appreciate any help resolving this issue. I am using Mathematica version 10.4.0.0.

• Sure seems like buggy behavior to me. I see similar results in Mathematica 6, so it's not anything new, at least... – Simon Rochester Nov 17 '16 at 6:39
• Maybe this Mathoverflow question and answer are useful, especially if you're just interested in the values and not necessarily have to use Mathematica. – Jules Lamers Nov 18 '16 at 13:08
• Thank you very much for your suggestions! I tried to use Maple for my calculation and I got the result that I was looking for. – Saya Nov 19 '16 at 8:54

Yes, it seems to be a bug, and partly related to the bug described in SiegelTheta throws errors from calling Range with complex arguments. For reference, the old implementation by Keiper in earlier versions (I used 5.2 for this example) is able to deal with this:

N[SiegelTheta[{{2 I, I, -3, 3}, {I, 2 I, 0, -3}, {-3, 0, 2 I, I}, {3, -3, I, 2 I}},
{0, 0, 0, 0}], 120]
1.0225349362774725133736050639011899697771634930261070909537271092899151358919942752331766696357309426506308161902561241867315171877661039299108073325227634120.15051499783199 + 0120.1408368426457*I


IIRC, the SiegelTheta[] implementation uses the Deconinck-van Hoeij method internally, so I'm not sure why it's slipping up this way.

Edit: After the comment by @Saya I can only say that this looks like a bug to me. I am afraid my reply does not answer anything, and only affirms that there seems to be a problem. In any case, here is my (slightly adapted) original reply, followed by an addendum.

The following is too long to fit in a comment, and might be useful, even though I'm not sure it qualifies as an answer. In Mathematica 11.0, when I quit the kernel, and then evaluate

P = {{2 I, I, -3, 3}, {I, 2 I, 0, -3}, {-3, 0, 2 I, I}, {3, -3, I, 2 I}};
N[SiegelTheta[P, {0, 0, 0, 0}], 50]
N[SiegelTheta[P, {0, 0, 0, 0}], 60]
N[SiegelTheta[P, {0, 0, 0, 0}], 70]
N[SiegelTheta[P, {0, 0, 0, 0}], 100]


I get the same strange sequence of outputs that you get, together with a bunch of errors of the form

Range::range: Range specification in Range[-I,I] does not have appropriate bounds.


for various Ranges, always including an I somewhere. I have no idea where these Ranges with imaginary arguments appear in the evaluation. However, when evaluating the same lines a second time the error messages disappear and Mathematica seems to give sensible results:

(* Out *)
4.8958265851687817320116929843990980640568162474355*10^46 + 0.*10^-4 I
4.89582658516878173201169298439909806405681624743553153865985*10^46 + 0.*10^-14 I
4.895826585168781732011692984399098064056816247435531538659849186987472*10^46 + 0.*10^-24 I
4.895826585168781732011692984399098064056816247435531538659849186987471892518027397637334165655318210*10^46 + 0.*10^-54 I


(Just evaluating SiegelTheta[P, {0, 0, 0, 0}] gives back the full form, with P substituted, but no error. The problem has to lie in the way Mathematica computes the numeric value, but I'm affraid I can't offer any more insight.)

Edit. Based on your comment let's go a bit further and next evaluate

N[SiegelTheta[P, {0, 0, 0, 0}], 50]
N[SiegelTheta[P, {0, 0, 0, 0}], 60]
N[SiegelTheta[P, {0, 0, 0, 0}], 70]
N[SiegelTheta[P, {0, 0, 0, 0}], 100]
N[SiegelTheta[P, {0, 0, 0, 0}], 120]


Up to 100 digits, which is what we asked for before, Mathematica just returns the appropriate number of digits of the most precise value that it computed so far. The last line, however, once more yields Range errors, and the final value is different again:

(* Out *)
4.8958265851687817320116929843990980640568162474355*10^46 + 0.*10^-4 I
4.89582658516878173201169298439909806405681624743553153865985*10^46 + 0.*10^-14 I
4.895826585168781732011692984399098064056816247435531538659849186987472*10^46 + 0.*10^-24 I
4.895826585168781732011692984399098064056816247435531538659849186987471892518027397637334165655318210*10^46 + 0.*10^-54 I

(* Let me omit the Range::range's *)

4.02565781844434120871609151731257530059095168812684514304059622770541639174945549458388772880374954035010464735943363847*10^57 + 0.*10^-63 I


This pattern continues; it always returns the most precisely computed value (to the numbers of digits asked for), and when increasing the precision we get the Range` errors together with a new value that is different (often by many orders of magnitude) from the most precise value found before...

• When evaluating the same lines a second time I also get the result evaluated with the highest precision, N[#,100], whichever of N[#,50], N[#,60],N[#,70] and N[#,100] I try to evaluate. I speculate that Mathematica may temporarily retain the result with the highest precision and reuse it in calculating the same in order to increase speed of calculation, so we may get wrong results in evaluating a second time. – Saya Nov 17 '16 at 14:09
• @Saya you're right; I have adjusted my answer accordingly. I'm afraid I can't help! – Jules Lamers Nov 17 '16 at 14:51