Guest User

Untitled

a guest
Sep 24th, 2016
38
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include <stdio.h>
  2. #include <math.h>
  3. #include <stdlib.h>
  4.  
  5. #define N_1 10
  6.  
  7. double chebyshev(int n, double x){
  8.     switch(n){
  9.         case 0:
  10.             return 1.0;
  11.             break;
  12.         case 1:
  13.             return x;
  14.             break;
  15.         default:
  16.             return 2.0*x*chebyshev(n-1, x)-chebyshev(n-2, x);
  17.             break;
  18.     }
  19. }
  20.  
  21. double shifted_chebyshev(int n, double x){
  22.     return chebyshev(n, 2.0*x-1.0);
  23. }
  24.  
  25. double interval1(double x){
  26.     static const double pj_1[N_1] = {
  27.                             -4.1658081333604994241879E11,
  28.                             +1.2177698136199594677580E10,
  29.                             -2.5301823984599019348858E10,
  30.                             +3.1984354235237738511048E8,
  31.                             -3.5377809694431133484800E8,
  32.                             -3.1398660864247265862050E5,
  33.                             -1.4299841572091610380064E6,
  34.                             -1.4287072500197005777376E4,
  35.                             -1.2831220659262000678155E3,
  36.                             -1.2963702602474830028590E1
  37.     };
  38.     static const double qj_1[N_1] = {
  39.                             -1.7934749837151009723371E11,
  40.                             +9.8900934262481749439886E10,
  41.                             -2.8986272696554495342658E10,
  42.                             +5.4229617984472955011862E9,
  43.                             -7.0108568774215954065376E8,
  44.                             +6.4698830956576428587653E7,
  45.                             -4.2648434812177161405483E6,
  46.                             +1.9418469440759880361415E5,
  47.                             -5.5648470543369082846819E3,
  48.                             +7.6886718750000000000000E1
  49.     };
  50.     static const double x0 = 0.372507410781366634461991866580;
  51.  
  52.     int j;
  53.     double result, numerator, denominator, sum;
  54.  
  55.     result = log(x/x0);
  56.     numerator = 0.0;
  57.     denominator = 0.0;
  58.     sum = 0.0;
  59.  
  60.     for(j=0; j<N_1; j++){
  61.         numerator += pj_1[j] * shifted_chebyshev(j, x/6.0);
  62.         denominator += qj_1[j] * shifted_chebyshev(j, x/6.0);
  63.     }
  64.     sum = numerator / denominator;
  65.     sum *= (x-x0);
  66.     result += sum;
  67.  
  68.     return result;
  69. }
  70.  
  71. double ei(double x){
  72.     if( x < 0.0 ){
  73.         printf("Argument of ei(x) must be positive\n");
  74.         exit(1);
  75.     }
  76.  
  77.     if( x <= 6.0 ){
  78.         return interval1(x);
  79.     }
  80.  
  81.     printf("Argument range not yet implemented for ei(x).\n");
  82.     exit(1);
  83. }
  84.  
  85. int main(void){
  86.     double x = 5.0;
  87.     printf("ei(%g)=%g\n", x, ei(x));
  88.    
  89.     return 0;
  90. }
RAW Paste Data