Fine tuning compiled code that computes dilogarithm function



As an exercise in writing a good Compile function, I want to do the simple task of coding a routine that outputs the real part of the dilogarithm function reLi2[z], given a complex number z as input. I would like some feedback in my code below that does the job.

I am following the strategy outlined in Celestial Mech. Dynam. Astron. 62 (1): 93–98. The idea is to divide the complex $z$-plane into four regions, as I show in the figure below, and to implement a different approximation in each region.

For simplicity on this site, I only carry out the task of evaluating only on the real $z$-axis.

complex plane regions

Edit: Based on ybeltukov's comment, I have updated the code (in ways I'd like to understand better) to make it so that SetSystemOptions["CompileOptions" -> "CompileReportExternal" -> True]; doesn't complain.

  1. In region 1, approximate the dilogarithm by the defining sum: $\operatorname{Li}_2(z) \approx \sum_{k=1}^\infty x^k/k^2$. Here is the code:

    realRegion1 = Compile[{{x, _Real}}, Sum[x^k/k^2, {k, 1., 23}]];
  2. In region 2, approximate the dilogarithm by evaluating the integral $\operatorname{Li}_2(z) \approx -\int_0^1 \frac{\ln(1-z t)}{t}\,dt$ by Gaussian quadrature (9 divisions):

    With[{div = 9},
      With[{x = Sort[formalX /. Solve[LegendreP[div, formalX] == 0, formalX] // N]},
          y = Chop[x/2 + 1/2], 
          w = Table[2/((1 - x[[i]]^2)*Derivative[0, 1][LegendreP][div, x[[i]]]^2), 
                {i, 1, div}]},
            realRegion2expr =
                      Chop[w[[i]]] Log[1 - 2 y[[i]] var + y[[i]]^2 var^2]/y[[i]], 
                      {i, 1, div}]},
            realRegion2 = Compile[{{var, _Real}}, realRegion2expr]]
  3. In region 3, apply the dilogarithm identity $\operatorname{Li}_2(z) = -\underbrace{\operatorname{Li}_2(1-z)}_\text{region I} - \ln(z)\ln(1-z)+\pi^2/6$, where the dilogarithm in the RHS is to be evaluated in region I.

    realRegion3 = Compile[{{x, _Real}}, 
     If[x == 1, Pi^2/6., -realRegion1[1 - x] - Log[x] Log[1 - x] + Pi^2/6.],
       CompilationOptions -> {"InlineExternalDefinitions" -> True}];
  4. In region 4, apply the dilogarithm identity $\operatorname{Li}_2(z) = -\underbrace{\operatorname{Li}_2(1/z)}_\text{region I,II,III} - \frac{1}{2}\ln^2(-z)-\pi^2/6$, where the dilogarithm in the RHS is to be appropriately evaluated in region I, II or III depending on the value of $1/z$.

    So first, I need to put together the functions realRegion1, realRegion2 and realRegion3 so that it correctly evaluates on the real line segment $-1 \leq z \leq +1$ appropriately:

    realSegment = 
      Compile[{{x, _Real}},
        If[-0.5 <= x <= 0.5, realRegion1[x],
          If[x <= 0, realRegion2[x], realRegion3[x]]
        ], CompilationOptions -> {"InlineExternalDefinitions" -> True}

And now, I can do region IV (and also including everywhere else) for my final compiled function

reLi2 = 
  Compile[{{x, _Real}}, 
    Re@If[-1. <= x <= 1., 
        -realSegment[1/x] - Pi^2/6 - 1/2*(1/4 Log[x^2]^2 - Arg[-x]^2)], 
  CompilationOptions -> {"InlineExternalDefinitions" -> True}]

The output is quite satisfactory. You can compare Plot[reLi2[x], {x, -5, 5}] with Plot[Re @ PolyLog[2, x], {x, -5, 5}].

However, I have no idea if I have compiled my function correctly for substantial increase in speed. I appreciate any feedback, no matter how minor, on my code.


Posted 2014-10-11T12:39:01.570

Reputation: 18 597

You can increase the speed in the first region with Sum[x^k/k^2, {k, 1., 23}]. You can also use options RuntimeAttributes -> {Listable}, Parallelization -> True if you want to calculate the function for a big list of arguments. – ybeltukov – 2014-10-11T13:32:30.237

1Actually you don't use compile benefits since SetSystemOptions["CompileOptions" -> "CompileReportExternal" -> True] reports that your main functions cannot be compiled and will be evaluated externally. Please try to rewrite your code step by step to avoid uncompiled evaluations. – ybeltukov – 2014-10-11T13:43:10.593

Something strange is going on. When use SetSystemOptions[...] as you suggested, I get an error that region1, region3, etc.. all can't be compiled. But this is not true. For example, I can compile region1 just fine by replacing Function with Compile... strange.. – QuantumDot – 2014-10-11T14:09:21.470

Ok, I will fix this... – QuantumDot – 2014-10-11T14:12:10.567

I finished fixing the code so that SetSystemOptions[...] no longer complains. Many thanks to @ybeltukov for pointing this out, but I'm not entirely sure how it will benefit me. Any other suggestions would be most helpful. – QuantumDot – 2014-10-11T22:32:38.533

1@ybeltukov You know something interesting? My original code Sum[x^k/k^2., {k, 1, 23}] went faster than your suggestion Sum[x^k/k^2, {k, 1., 23}]. And Sum[x^k/k^2, {k, 1, 23}] goes even faster still! Maybe raising numbers to power exact 2 is better than raising to approximate 2, and performing sum with exact numbers is faster than with approximate numbers... – QuantumDot – 2014-10-11T22:41:34.037

Good work! It is really interesting... For further speed up you can add RuntimeOptions -> "Speed", CompilationTarget -> "C". After it my test shows, that your function is 10 faster then built-in Re@PolyLog[2, x]. Is it good enough? – ybeltukov – 2014-10-11T22:47:16.627

2could you post an answer to your question if you believe you now have a fast solution? – chris – 2015-06-07T20:19:42.410


I have a feeling you'll want to see this paper.

– J. M.'s ennui – 2015-06-07T20:29:48.960



Here is a compiled routine for evaluating the dilogarithm $\operatorname{Li}_2(x)$ for real $x$, using the fourth-order series in Morris's paper (linked in the comments), along with the use of functional equations to bring the argument into a range where the series can be efficiently evaluated. (It is due to these functional equations that the code is a bit on the long side.)

dilog = With[{eps = $MachineEpsilon, pi26 = N[π^2/6]},
             Compile[{{x, _Real}},
                     Module[{c = 1., s = 3., j = 1, k = 1, l0 = 1., l1 = 1., l, t, xx},
                            If[x == 1., Return[pi26]];
                            If[x > 0.5, l0 = Log[x]]; If[x < 1., l1 = Internal`Log1p[-x]];
                            xx = Which[x < -1., 1./(1. - x),
                                       -1. <= x < 0., x/(x - 1.),
                                       0.5 < x < 1., 1. - x,
                                       1. < x < 2., 1. - 1./x,
                                       2. <= x, 1./x,
                                       True, x];
                            While[l = k + 2 j + 1; c *= xx;
                                  t = c/(k l); k = l; s += t; j++;
                                  Abs[t] >= eps Abs[s]];
                            s = xx/(xx + 1) s - 2 (xx - 1)/(xx + 1) Internal`Log1p[-xx];
                            Which[x < -1., s - pi26 + l1 (0.5 l1 - Log[-x]),
                                  -1. <= x < 0., -0.5 l1^2 - s,
                                  0.5 < x < 1., pi26 - s - l0 l1,
                                  1. < x < 2., pi26 + l0 (0.5 l0 - Log[x - 1.]) + s,
                                  2. <= x, 2. pi26 - 0.5 l0^2 - s,
                                  True, s]],
                     RuntimeAttributes -> {Listable}]];

One might notice the use of the function Internal`Log1p[] for evaluating $\log(1+x)$, which only became recently available. For versions that do not yet have this function, a trick like the one in this answer could be used.

As is conventional for $x>1$, only the real part is returned; You can add the imaginary part $-i\pi\log x$ yourself if needed.

The routine dilog does well in accuracy:

Plot[{Re[PolyLog[2, x]], dilog[x]}, {x, -30, 30},
     PlotLegends -> {PolyLog[2, x], "dilog(x)"}, 
     PlotStyle -> {AbsoluteThickness[6], AbsoluteThickness[2]}]

PolyLog vs. dilog

Plot[-Log10[Abs[dilog[x]/Re[PolyLog[2, x]] - 1]], {x, -30, 30}, 
     PlotLabel -> "Relative Error", PlotRange -> All]

relative error of dilog

(The dip to the right in the plot of the relative error corresponds to the zero of $\Re\operatorname{Li}_2(x)$ at $x_\ast\approx 12.59517$. You can compute $x_\ast$ with FindRoot[π^2/3 - Log[z]^2/2 - PolyLog[2, 1/z], {z, 12}].)

For a speed comparison:

zs = N[Range[-10, 10, 1*^-4]];
AbsoluteTiming[Re[PolyLog[2, zs]];]
   {4.01328, Null}

   {0.150729, Null}

J. M.'s ennui

Posted 2014-10-11T12:39:01.570

Reputation: 115 520