# Code to calculate the partial theta function

The version of the partial theta function I want to compute is

$O(z) = \sum\limits^\infty_{n=0}\exp(-(z+n\pi)^2)$

Does anyone have or know of code to efficiently compute $O(z)$? There are two papers dealing with this computation, but can't find any code.

Hiary, G. (2008). A nearly-optimal method to compute the truncated theta fcn, its derivatives and integrals. Annals of Mathematics, Second Series, Vol 174, No. 2, pp 859-889.

Kuznetsov, A (2015). Computing the truncated theta function via Mordel-Integral. Math comp 84, pp 2911-2926

• I'm guessing this is not a mathematica question, but if it is start by whats wrong with NSum[ Exp[-(z + n Pi)^2 ], {n, 0, Infinity}] – george2079 Jan 23 '18 at 17:08
• Is there an I missing? – Michael E2 Jan 23 '18 at 18:03

θ[z_?NumericQ] := NSum[Exp[-(z + n Pi)^2], {n, 0, Infinity}]


Using a partial sum

θp[z_] := Sum[Exp[-(z + n Pi)^2], {n, 0, 100}]


Comparing the plots

Plot[{θ[z], θp[z]}, {z, -2.5, 2.5},
PlotStyle -> {{AbsoluteThickness[4], Green}, Red},
PlotLegends -> Placed["Expressions", {.2, .75}]]


• just think someone got a journal paper out of that :-). (Something tells me there is more to this.) – george2079 Jan 23 '18 at 17:42

The sum $O(z) = \sum\limits^\infty_{n=0}\exp(-(z+n\pi)^2)$ is extremely rapidly convergent, unlike the sums in the references cited in the question.

Pick an accuracy eps. Then one can sum to that accuracy as follows:

(* solve for limits of summation *)
n /. Simplify[Solve[Exp[-(z + n Pi)^2] == eps, n, Reals], 0 < eps < 1]
(*  {-((z + Sqrt[-Log[eps]])/π), (-z + Sqrt[-Log[eps]])/π}  *)

ClearAll[sum];
sum[z_?NumericQ, accgoal_: Automatic] := Module[{n, eps},
eps = 10^-(accgoal /. Automatic -> $MachinePrecision); n = Range @@ (* range of terms to sum *) N@Clip[ Floor@{-((z + Sqrt[-Log[eps]])/π), (-z + Sqrt[-Log[eps]])/π} + {0, 1}, {0, Infinity}]; Total[Exp[-(z + n Pi)^2]] (* vectorized sum for speed *) ]; Plot[sum[z], {z, -30, 10}]  The number of terms in the sum for machine precision accuracy is about 4 or 5: 2*Sqrt[-Log[$MachineEpsilon]]/Pi
(*3.82203  *)