C#で対数積分関数li(x)を自前で計算する(x≧2)

$$\text{li}(x):=\lim_{\epsilon\to 0+}\left(\int_0^{1-\epsilon}\frac{dt}{\ln t}+\int_{1+\epsilon}^x\frac{dt}{\ln t}\right)\ \ \ \ \ (x>1)$$

を計算するプログラムを書いてみた。

上記で紹介されているラマヌジャンの式を使うと高速に収束するとのことなので、これを用いる。
$$\text{li}(x)=\gamma+\ln\ln x + \sqrt{x}\sum_{n=1}^\infty\frac{(-1)^{n-1}(\ln x)^n}{n!2^{n-1}}\sum_{k=0}^{\lfloor(n-1)/2\rfloor}\frac{1}{2k+1}$$

li.cs

using System;

class LiFunc
{
    static double gamma_constant = 0.57721566490153286061;

    static double Li(double x)
    {
        double lnX = Math.Log(x);
        double a = gamma_constant + Math.Log(lnX);
        double sum1=0;
        double prevSum1=0;
        double sum2=0;
        double factor = lnX;
        int prevK = -1;
        int endK;
        for(int n=1;n<=100000;n++)
        {
            if (n>1)
            {
                factor *= -lnX*0.5/n;
            }
            endK = (n-1)/2;
            for(int k=prevK+1;k<=endK;k++)
            {
                sum2 += 1.0/(2*k+1);
            }
            prevK = endK;

            sum1 += factor*sum2;
            if (prevSum1 == sum1)
            {
                // doubleの範囲で値が変化しなくなったときのnを出力(参考用)
                Console.Write("exit at n=");
                Console.WriteLine(n);
                break;
            }
            prevSum1 = sum1;
        }
        return a + sum1*Math.Sqrt(x);
    }

    [STAThread]
    static void Main(string[] args)
    {
        Console.WriteLine(Li(2));
        Console.WriteLine(Li(1.1));
        Console.WriteLine(Li(1.4513692348)); // 計算結果が 0 に近いため計算誤差が累積する。
        Console.WriteLine(Li(10000000000000));
    }
}

出力結果

exit at n=14
1.04516378011749
exit at n=10
-1.6757728303192
exit at n=12
-2.23837337554045E-10
exit at n=61
346065645810.05

python (mpmath)で検算

calc_li.py

import mpmath
mpmath.mp.dps = 17

for x in [
    2,
    1.1,
    1.4513692348,
    10000000000000,
    ]:
    print(mpmath.li(x))

出力結果

1.0451637801174928
-1.6757728303191988
-2.2383724226679356e-10
346065645810.05034

x=1.4513692348
以外の例はdoubleの桁数の範囲でmpmathの計算結果を丸めたものと一致しているようなので、$x\geq 2$のときはわりと実用的に使えるレベルかなと。