#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;
}