Tests of lerchphi()

Sergej V. Aksenov
Department of Microbiology and Immunology
University of Michigan
Ann Arbor, MI 48109, USA
http://aksenov.freeshell.org/

Ulrich D. Jentschura
Institute of Theoretical Physics
Dresden University of Technology
01062 Dresden, Germany
jentschura@physik.tu-dresden.de

Regular cases in output from program test.c (cases 5-12), saved to test.ref

First, some definitions...

form[x_]  :=  
      PaddedForm[x, {34, 31},  DigitBlock  ->  5,  
        NumberSeparator  ->  " "]

shortform[x_]  :=  
      PaddedForm[x, {22, 16},  DigitBlock  ->  5,  
        NumberSeparator  ->  " "]

Attributes[form] = {Listable};  
    Attributes[shortform] = {Listable};

Example cases 5 through 12:

We ask Mathematica to calculate in extended precision (quad-double):

testcases = {{0.99999`64, 2.0`64, 1.0`64}, {-0.99999`64, 2.0`64, 
            1.0`64}, {0.99999`64, 2.0`64, 1.0`64 10^3}, {0.3`64, 
            2.0`64, -4.5`64}, {0.00001`64, 2.0`64, 
            1.0`64}, {-0.000063`64, 2.0`64, 
            1.0`64}, {3.4709929976435479`64 10^(-6), 1.0`64, 
            1.5172413793103448`64}, {0.0003`64, 
            2, -3.00000000000001`64}};

The table of the input parameters z, s, v reads:

Join[{{"z", "s", "v"}}, shortform[testcases]] // TableForm

$MinPrecision = 64

$Version

Timing[
        testres = Map[Apply[LerchPhi, #] &, testcases]][[1]]

A typical execution time for the C++ implementation (file lerchphimp.cc) is 0.6 seconds on a Pentium II. This implementation uses optimized arithmetic with 64 significant decimals ("quad-double" precision) (go to web page). Note that the computational overhead in Mathematica caused by handing the arguments over to the routine LerchPhi should be small. Our LerchPhiCNCT implementation in Mathematica is not applicable to the parameter region z < 0, and thus cannot be applied to the example calculations 6 and 10 (see below).

factors = {1.`64, 10.`64, 10000.`64, 10.`64, 1.`64, 10.`64, 10.`64, 
          1.0`64 10^(-17)};

casenum = Table[i, {i, 5, 12}];

ResTable  =  
      TableForm[Transpose[{casenum, form[factors testres]}]]

As shown below, cases 7 and 12 are interesting (the above results are correct).

Export["test.ref", ResTable, "Table"];

The Problematic Cases 7 and 12

Case 7: The result from the C program (lerchphi.c) in standard double precision x 10^4 is     9.59714 89709 9796    which is - on the level of one part in 10^12 - different from the above result of  9.59714 89709 9654...

Case 12: The result from the C program (lerchphi.c) in standard double precision x 10^(-17) is 2.5880201290103731 which is in disagreement with the result of 2.70000000000000000111 in the table above.

These are the only two cases where the standard-precision routine (lerchphi.c) written in C (double, 16 digits) is not sufficient to calculate the result to the desired precision (14 digits).
 
 An explanation for the numerical problems occurring in both of the above cases is given in the program documentation file  (lphidoc.pdf).
 
The discrepancy disappears if the algorithm is augmented by extended-precision arithmetic either in Mathematica (both with the built-in routine LerchPhi and also with our implementation LerchPhiCNCT), and it also disappears when using the multiprecision C++ routine (lerchphimp.cc).

Load CNCT algorithm from our Mathematica package.

<<  ../Mathematica/lerchphi.m

Require more accuracy in the result.

SetOptions[LerchPhiCNCT, CNCTRelativeAccuracy -> 32]

In a C program one can use a multiprecision software library, for example by Bailey et al. (go to web page) or by Briggs (go to web page).

Case 7

The result from the C program (lerchphi.c) in standard precision x 10^4 is     9.59714 89709 9796.

The result of Mathematica's built-in LerchPhi is:

form[10^4 testres[[7 - 5 + 1]]]

The result obtained with our implementation LerchPhiCNCT (in Mathematica) coincides with the built-in LerchPhi routine:

form[10^4 Apply[LerchPhiCNCT, testcases[[3]]]]

The same is true for the output of lerchphimp.cc (result not shown).

Case 12

The result from the C program in standard precision x 10^(-17) is     2.58802 01290 1037.

The result of Mathematica's built-in LerchPhi is:

form[testres[[12 - 5 + 1]]]

The result obtained with our implementation LerchPhiCNCT (in Mathematica) coincides with the built-in LerchPhi routine:

form[10^4 Apply[LerchPhiCNCT, testcases[[8]]]]

The same is true for the output of lerchphimp.cc (result not shown).