      function ei=ei(x)
      f=-1.e-99;
      if(abs(x-10) < 12)
        j=fix(10+2*abs(x)) ;
        f=(j+1)^(-2)       ;
        while (j > 0)
         f=(f*j*x+1)*j^(-2);
         j=j-1             ;
        end
        ei=f*x+log(1.781072418*abs(x));
        return
      else
        j=fix(5+20/abs(x)) ;
        f=x                ;
        while (j > 0)
         f=1/(1/f-1/j)+x   ;
         j=j-1             ;
        end
        ei=exp(x)/f        ;
        return
      end
%  from An Atlas of Functions, Spanier and Oldham
