# Number of divisors visualized with the QPochhammer function, how to improve performance of code?

I have this code that is originally Jeffrey Stopple's code for the Riemann zeta function in the complex plane. Because I discovered yesterday that the number of divisors can be generated with the $q$-Pochhammer symbol (QPochhammer), and since Mathematica shows a plot of QPochhammer, I thought that plotting it would be fun.

Here is the code that needs improvement:

Show[Graphics[RasterArray[
Table[Hue[Mod[3 Pi/2 +
Arg[Sum[(s + I t)^(n - 1)*(QPochhammer[(s + I t)^(n + 1), (s + I t)]/
QPochhammer[(s + I t)^(n), (s + I t)]), {n, 1, 100}]],
2 Pi]/(2 Pi)], {t, -1.1, 1.1, .05}, {s, -1.1, 1.1, .05}]]],
AspectRatio -> Automatic]


And the code for the number of divisors:

CoefficientList[
Series[Sum[
x^(n - 1)*(QPochhammer[x^(n + 1), x]/QPochhammer[x^(n), x]), {n, 1,
104}], {x, 0, 103}], x] (*_ Mats Granvik_,Jan 03 2015*)

1, 2, 2, 3, 2, 4, 2, 4, 3, 4, 2, 6, 2, 4, 4, 5, 2, 6, 2, 6, 4, 4, 2, 8, 3, 4, 4,
6, 2, 8, 2, 6, 4, 4, 4, 9, 2, 4, 4, 8, 2, 8, 2, 6, 6, 4, 2, 10, 3, 6, 4, 6, 2, 8,
4, 8, 4, 4, 2, 12, 2, 4, 6, 7, 4, 8, 2, 6, 4, 8, 2,...


The plot from the first program that I would like to improve:

Already if someone could post a plot with higher resolution, I would be glad. My computer is rather old.

## 2 Answers

It looks like you want to plot the phase-only information of a complex function. Using the following helper functions for plotting the phase-only information complex functions:

hue = Compile[{{z, _Complex}}, {Mod[3 π/2 + Arg[z],
2 π]/(2 π), 1, If[Abs[z] > 10^-3, 1, 0]},
CompilationTarget -> "C", RuntimeAttributes -> {Listable}];
ComplexPlotC[f_, {x0_, x1_, δx_}, {y0_, y1_, δy_}] :=
Image[hue[
f[Outer[Complex, Range[x0, x1, δx],
Range[y1, y0, -δy]]]]\[Transpose], ColorSpace -> Hue,
Magnification -> 1];
CCompileC[expr_] :=
Compile[{{z, _Complex}}, Evaluate[expr], CompilationTarget -> "C",
RuntimeAttributes -> {Listable}];


You can then compile your function and plot it (I'll only sum 10 terms, rather than 100, as using 100 would take quite a long time):

func = CCompileC[
If[Abs[z] >= 0.999, 0,
Sum[z^(n - 1)*(QPochhammer[z^(n + 1), z]/
QPochhammer[z^(n), z]), {n, 1, 10}]]];
ComplexPlotC[func, {-1 + 10^-6 RandomReal[], 1,
0.003}, {-1 + 10^-6 RandomReal[], 1, 0.003}]


which gives the following:

Here is a 100-term sum picture (open in separate tab to see slightly larger picture):

• damned: you scooped my 'brute force' solution by 2 minutes :-) – chris Jan 4 '15 at 21:20
• ah… but you stop at n=10? – chris Jan 4 '15 at 21:24
• @chris: Oops, good point, I changed it to 10 and totally forgot about it! Will edit. Ugh, now it's going to be really slow... – DumpsterDoofus Jan 4 '15 at 21:25
• well that's why it was taking sooo long on 40 cores :-) – chris Jan 4 '15 at 21:25
• Your plot is really nice! I like it. Thanks. – Mats Granvik Jan 4 '15 at 21:26

As this is a question, I feel justified in using a bit of heavy artillery. Here goes nothing...

In effect, what the OP seems to want to do is to evaluate

$$\sum_{n=1}^\infty \frac{(q^{n+1};q)_\infty}{(q^n;q)_\infty} q^{n-1}$$

(where $(a;q)_n$ is the $q$-Pochhammer symbol) by approximating it with its partial sums. However, there is a more streamlined expression for this function.

Using, for instance, formula 12 in the MathWorld page I linked to, we have

\begin{align*}\sum_{n=1}^\infty \frac{(q^{n+1};q)_\infty}{(q^n;q)_\infty} q^{n-1}&=\frac{(q^2;q)_\infty}{(q;q)_\infty}\sum_{n=0}^\infty \frac{\frac{(q;q)_\infty}{(q^{n+1};q)_\infty}}{\frac{(q^2;q)_\infty}{(q^{n+2};q)_\infty}} q^n\\&=\frac1{(q;q)_1}\sum_{n=0}^\infty \frac{(q;q)_n}{(q^2;q)_n} q^n\\&=\frac1{1-q}\sum_{n=0}^\infty \frac{(q;q)_n (q;q)_n}{(q^2;q)_n}\frac{q^n}{(q;q)_n}\end{align*}

and the infinite sum can now be recognized as the $q$-hypergeometric function $$\frac1{1-q}\,{}_2\phi_1\left({{q,q}\atop{q^2}}\mid q,q\right)$$

or in Mathematica syntax, QHypergeometricPFQ[{q, q}, {q^2}, q, q]/(1 - q).

Another expression is derived by noting that

$$\frac{(q;q)_n}{(q^2;q)_n}=\frac{1-q}{1-q^{n+1}}$$

Using this identity instead in the derivation above yields

\begin{align*}\frac1{1-q}\sum_{n=0}^\infty \frac{(q;q)_n}{(q^2;q)_n} q^n&=\frac1{q}\sum_{n=1}^\infty \frac{q^n}{1-q^n}\end{align*}

i.e., a Lambert series, which is expressible in terms of the $q$-digamma function (compare the definitions)

$$\frac{\psi_q(1)+\log(1-q)}{q\log q}$$

or (QPolyGamma[1, q] + Log[1 - q])/(q Log[q]) in Mathematica syntax.

We now have two alternative expressions to choose from. The $q$-digamma function is faster to evaluate (at least in the tests I did), but the $q$-hypergeometric expression does not have the removable singularity at the origin.

If we are to use the $q$-digamma function expression, then, this precludes us from an Image[]-based approach such as the one in the OP and the other answer. We could use DensityPlot[] (since it does not visibly complain about singularities) with an appropriate RegionFunction setting, but a better way would be to use another classical technique, this time from Schwarz and Christoffel.

To be precise, what I'll do here is to use the conformal map from the unit square to the unit disk within DensityPlot[] and then perform an ImageForwardTransformation[] (like in the one done here) to finally obtain the image on the unit disk. As mentioned in that other answer, the appropriate mapping function is JacobiSN[z, 1/2] JacobiDC[z, 1/2]. The built-in Jacobi elliptic functions are a bit slow for image transformation purposes, so I have written a compiled function that evaluates this mapping for complex arguments:

sdc = With[{c1 = 1./2678, c2 = 1101., c3 = 110.2, c4 = 5677./60,
c5 = 1577., c6 = 1479., c7 = 21739./60, d1 = 1./338, d2 = 51.,
d3 = 1./15, d4 = 1709., d5 = 262.25, d6 = 287., d7 = 0.2,
d8 = 1591./3, d9 = 102.75},
Compile[{{z, _Complex}}, Module[{dc, k, s, zs, zs2},
k = If[z == 0., 0, Max[0, Ceiling[Log2[4. Abs[z]]]]];
zs = z/(2^k); zs2 = (zs/2)^2;
{s, dc} = {(zs (1. - c1 zs2 (c2 + zs2 (c3 - c4 zs2))))/
(1. + c1 zs2 (c5 - zs2 (c6 - c7 zs2))),
(1. + d1 zs2 (d2 + d3 zs2 (d4 + d5 zs2)))/
(1. - d1 zs2 (d6 + d7 zs2 (d8 + d9 zs2)))};
Do[{s, dc} = {(2. s dc)/(1. + s^2 dc^2),
(dc^2 - 0.5 s^2)/(1. - s^2 dc^2)}, {k}];
s dc], RuntimeAttributes -> {Listable}]]


We now generate the image in two stages. First, here's the "raw" DensityPlot[]:

With[{ω = N[3 Beta[5/4, 5/4]/2], ε = 1/100},
sqrimg = DensityPlot[With[{q = sdc[ω (x + I y)]},
Arg[(QPolyGamma[1, q] + Log[1 - q])/(q Log[q])]],
{x, -1 + ε, 1 - ε}, {y, -1 + ε, 1 - ε},
ColorFunction -> (Hue[Mod[#/(2 π) + 3/4, 1]] &),
ColorFunctionScaling -> False, Frame -> None,
ImageSize -> Large, MaxRecursion -> 1,
PlotPoints -> 95, PlotRangePadding -> None]]


Some features should already be recognizable in this distorted image. Now, undo the Schwarz-Christoffel distortion like so:

With[{ω = N[3 Beta[5/4, 5/4]/2]},
ImageForwardTransformation[sqrimg, With[{z = ω (#[[1]] + I #[[2]])},
Through[{Re, Im}[sdc[z]]]] &, Background -> 1.,
DataRange -> {{-1, 1}, {-1, 1}}]]


Beautiful. If need be, you can always increase the setting for PlotPoints and MaxRecursion, but the picture already looks reasonably good to me. (Note that the jittery portion near $q=1$ in the other images posted is completely absent.)

• (A similar technique was used to render my current Gravatar, with the Rogers-Ramanujan continued fraction as the function being plotted.) – J. M. will be back soon Aug 1 '15 at 10:50
• These functions are not so familiar to me although I used them in my question, but your answer appears to be pretty rigourous. – Mats Granvik Aug 1 '15 at 10:52
• I try my best when it comes to dealing with special functions. :) – J. M. will be back soon Aug 1 '15 at 10:53