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

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]},
With[{
y = Chop[x/2 + 1/2],
w = Table[2/((1 - x[[i]]^2)*Derivative[0, 1][LegendreP][div, x[[i]]]^2),
{i, 1, div}]},
With[{
realRegion2expr =
-1/4.*Sum[
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[x],
-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.

• 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. Commented Oct 11, 2014 at 13:32
• Actually 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. Commented Oct 11, 2014 at 13:43
• @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... Commented Oct 11, 2014 at 22:41
• could you post an answer to your question if you believe you now have a fast solution? Commented Jun 7, 2015 at 20:19
• I have a feeling you'll want to see this paper. Commented Jun 7, 2015 at 20:29

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 = InternalLog1p[-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) InternalLog1p[-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 InternalLog1p[] 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]}]  Plot[-Log10[Abs[dilog[x]/Re[PolyLog[2, x]] - 1]], {x, -30, 30}, PlotLabel -> "Relative Error", PlotRange -> All]  (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}

AbsoluteTiming[dilog[zs];]
{0.150729, Null}
`