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

9

4

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:

Divisors from QPochhammer

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

Mats Granvik

Posted 2015-01-04T19:03:15.737

Reputation: 1 007

Answers

12

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:

enter image description here

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

enter image description here

DumpsterDoofus

Posted 2015-01-04T19:03:15.737

Reputation: 11 447

damned: you scooped my 'brute force' solution by 2 minutes :-) – chris – 2015-01-04T21:20:17.423

ah… but you stop at n=10? – chris – 2015-01-04T21:24:07.273

@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 – 2015-01-04T21:25:04.470

well that's why it was taking sooo long on 40 cores :-) – chris – 2015-01-04T21:25:33.720

Your plot is really nice! I like it. Thanks. – Mats Granvik – 2015-01-04T21:26:30.103

@MatsGranvik: It's also slightly wrong, as the sum only goes up to 10, rather than 100. – DumpsterDoofus – 2015-01-04T21:27:38.353

@DumpsterDoofus It does not matter. Your coding skills are amazing. I like the image anyways. – Mats Granvik – 2015-01-04T21:37:14.983

@MatsGranvik: The helper functions I used are actually something I wrote a long time ago, and only needed a slight modification. I'll post the 40-term sum picture when I get the parallelization to work correctly. – DumpsterDoofus – 2015-01-04T21:41:18.397

its got more funny loops... – chris – 2015-01-04T22:01:22.287

@chris: A pic with 2.5 times as many loops will hopefully finish up tonight... – DumpsterDoofus – 2015-01-04T22:33:24.383

Since QPochhammer is not compilable, you're probably better off not to compile func. hue, of course, can still freely be compiled. – Oleksandr R. – 2015-01-04T23:11:12.900

@OleksandrR.: Weird, it never gave any "proceeding with uncompiled evaluation" warnings when I tested the function inside the unit circle. It did give "proceeding with uncompiled evaluation" warnings when $|z|>1$, which is why I added the If condition. – DumpsterDoofus – 2015-01-04T23:26:03.033

10

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

distorted density plot of the Lambert series

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}}]]

The final image. :D

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

J. M.'s ennui

Posted 2015-01-04T19:03:15.737

Reputation: 115 520

(A similar technique was used to render my current Gravatar, with the Rogers-Ramanujan continued fraction as the function being plotted.) – J. M.'s ennui – 2015-08-01T10:50:48.143

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 – 2015-08-01T10:52:27.490

I try my best when it comes to dealing with special functions. :) – J. M.'s ennui – 2015-08-01T10:53:31.420