SiegelTheta throws errors from calling Range with complex arguments

Question

Bug introduced in 8.0 or earlier and persists through 12.1 or later. Fixed in 13.2.0 or earlier

---

This may or may not be related to the bug reported in this question.

I was trying to verify the results of this challenge over on Code Golf with the following code:

p = {{5., 2., 0., 0.}, {2., 5., 2., -2.}, {0., 2., 5., 0.}, {0., -2., 0., 5.}}
SiegelTheta[I*p/Pi, ConstantArray[0, Length@p]]

While I do get the correct result, SiegelTheta throws a bunch of of messages via Range along the lines of:

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

I've verified this with all of the example matrices in the linked challenge - the larger the matrix, the more of these messages are thrown. The same happens when I generate random positive definite matrices myself with the following code:

n = 5;
m = RandomChoice[{-1, 1}, {n, n}]
p = N@IdentityMatrix[n] + m\[Transpose].m;
PositiveDefiniteMatrixQ@p
SiegelTheta[I*p/Pi, ConstantArray[0, n]]

I have contacted Wolfram Support and they confirmed that this is a bug. If anyone can reproduce it in earlier Mathematica version, feel free to update the header.

Answer 1

Status update: in 13.2.0 we get

The code is provided below -see also OP:

p = {{5., 2., 0., 0.}, {2., 5., 2., -2.}, {0., 2., 5., 0.}, {0., -2., 0., 5.}}
SiegelTheta[I*p/Pi, ConstantArray[0, Length@p]]

n = 5;
m = RandomChoice[{-1, 1}, {n, n}]
p = N@IdentityMatrix[n] + m\[Transpose].m;
PositiveDefiniteMatrixQ@p
SiegelTheta[I*p/Pi, ConstantArray[0, n]]