11
$\begingroup$

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

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – ybeltukov Oct 11 '14 at 13:32
  • 1
    $\begingroup$ 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. $\endgroup$ – ybeltukov Oct 11 '14 at 13:43
  • 1
    $\begingroup$ @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... $\endgroup$ – QuantumDot Oct 11 '14 at 22:41
  • 2
    $\begingroup$ could you post an answer to your question if you believe you now have a fast solution? $\endgroup$ – chris Jun 7 '15 at 20:19
  • 1
    $\begingroup$ I have a feeling you'll want to see this paper. $\endgroup$ – J. M. is away Jun 7 '15 at 20:29
6
$\begingroup$

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}

AbsoluteTiming[dilog[zs];]
   {0.150729, Null}
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.