Numerical evaluation of the Exponential Integral Ei by rational Chebyshev approximations fails

Question

I am trying to evaluate the Exponential Integral $Ei(x)=-\int^{\infty}_{-x}\frac{e^{-t}}{t}dt$ for $x>0$ (interpreted as the Cauchy principal value) by using rational Chebyshev approximations, which can be found in a paper by Cody & Thacher "Chebyshev Approximations for the Exponential Integral $Ei(x)$". The paper is available online at ams.org. I am trying to use the equation on page 292 for the interval $0<x\le6$:

$$Ei(x)\simeq\log(x/x_0)+(x-x_0)\frac{\sum_{j=0}^{n}p_jT_j^{*}(x/6)}{\sum_{j=0}^{n}q_jT_j^{*}(x/6)}$$

Where $x_0\approx{}0.37$ is the zero of $Ei(x)$, $p_j$ and $q_j$ are the coefficients found in table II of the paper (I am using $n=9$) and $T_j^{*}(x)=T_j(2x-1)$ are the shifted Chebyshev polynomials. So I implemented the equation in a naive C program

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

#define N_1 10

double chebyshev(int n, double x){
    switch(n){
        case 0:
            return 1.0;
            break;
        case 1:
            return x;
            break;
        default:
            return 2.0*x*chebyshev(n-1, x)-chebyshev(n-2, x);
            break;
    }
}

double shifted_chebyshev(int n, double x){
    return chebyshev(n, 2.0*x-1.0);
}

double interval1(double x){
    static const double pj_1[N_1] = {
                            -4.1658081333604994241879E11,
                            +1.2177698136199594677580E10,
                            -2.5301823984599019348858E10,
                            +3.1984354235237738511048E8,
                            -3.5377809694431133484800E8,
                            -3.1398660864247265862050E5,
                            -1.4299841572091610380064E6,
                            -1.4287072500197005777376E4,
                            -1.2831220659262000678155E3,
                            -1.2963702602474830028590E1
    };
    static const double qj_1[N_1] = {
                            -1.7934749837151009723371E11,
                            +9.8900934262481749439886E10,
                            -2.8986272696554495342658E10,
                            +5.4229617984472955011862E9,
                            -7.0108568774215954065376E8,
                            +6.4698830956576428587653E7,
                            -4.2648434812177161405483E6,
                            +1.9418469440759880361415E5,
                            -5.5648470543369082846819E3,
                            +7.6886718750000000000000E1
    };
    static const double x0 = 0.372507410781366634461991866580;

    int j;
    double result, numerator, denominator, sum;

    result = log(x/x0);
    numerator = 0.0;
    denominator = 0.0;
    sum = 0.0;

    for(j=0; j<N_1; j++){
        numerator += pj_1[j] * shifted_chebyshev(j, x/6.0);
        denominator += qj_1[j] * shifted_chebyshev(j, x/6.0);
    }
    sum = numerator / denominator;
    sum *= (x-x0);
    result += sum;

    return result;
}

double ei(double x){
    if( x < 0.0 ){
        printf("Argument of ei(x) must be positive\n");
        exit(1);
    }

    if( x <= 6.0 ){
        return interval1(x);
    }

    printf("Argument range not yet implemented for ei(x).\n");
    exit(1);
}

int main(void){
    double x = 5.0;
    printf("ei(%g)=%g\n", x, ei(x));

    return 0;
}

But for the interval in question the results are completely off. For example I want to evaluate $Ei(5)\approx{}40.2$, but i am getting $\approx{}19.0654$. The paper claims a maximal relative error of $8\cdot10^{-19}$, so I guess I must be doing something wrong. I double checked the coefficients and my implementation of the equation but I don't see any mistake.

When I plot the relative error of the C calculated value compared to the real/expected value (as returned by the ExpIntegralEi function of Mathematica) I can see that only for the first interval ($0<x\le6$) the error is unreasonably high:

Plotted is $y(x)=\left|\frac{C\left(x\right)-M\left(x\right)}{M\left(x\right)}\right|$, where $C(x)$ is the value as returned by my C program and M(x) the same value as returned by Mathematica.

Any help would be much appreciated.

References

1. Cody, W. J., and Henry C. Thacher. "Chebyshev approximations for the exponential integral 𝐸𝑖 (𝑥)." Mathematics of Computation 23.106 (1969): 289-303.

Answer 1

There is a prime missing both in your code and in the expression in your question. In the original paper, the expression is: $$ \log(x/x_0) + (x-x_0) \frac{\sum^\prime_{0\leq j\leq n}p_j T_j^*(x/6)}{\sum^\prime_{0\leq j\leq n}q_j T_j^*(x/6)}. $$ This primed sum is very common and standard when dealing with Chebyshev series: it means that the first term of the sum, corresponding to the zeroth Chebyshev polynomial (the constant term), must be divided by two.

The origin of this modification is the identity ((25) on mathworld) $$ \sum_{1\leq k\leq m} T_i(x_k) T_j(x_k) = \frac{[i=j]m}{1 + [i=j\neq0]} $$ where $x_k$ are the roots of $T_m$.

After dividing $p_0$ and $q_0$ by $2$, I think it seems to work fine.

P.S. I would like to note that the standard, and really quite easy, way to evaluate Chebyshev series is with the Clenshaw recurrence, which is both much faster than your own implementation of Chebyshev polynomials (which takes time exponential in $n$ to evaluate $T_n$), and more accurate. In cases where coefficients have very large magnitudes, it can matter quite a bit, even though for some reason that doesn't seem to be the case here.