11

1

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.

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

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

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

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. – ybeltukov – 2014-10-11T13:32:30.2371Actually 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.593Something 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.470Ok, 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.5331@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.037Good 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.6272could you post an answer to your question if you believe you now have a fast solution? – chris – 2015-06-07T20:19:42.410

1

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

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