SiegelTheta provides a wrong value?

Question

Bug introduced in 6.0 or earlier and persists through 11.0.1 or later. Fixed in 13.2.0 or earlier.

---

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.

Answer 1

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.0225349362774725133736050639011899697771634930261070909537271092899151358919942752331766696357309426506308161902561241867315171877661039299108073325227634`120.15051499783199 + 0``120.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.

Answer 2

Just for completeness and clarity with 13.2.0 the results obtained are presented below

With

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

from the O.P we run

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]

and we obtain

It is evident that all the results are $O(1)$ as they should.

Answer 3

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...