forum.alglib.net

;===========================================
; E1 <= E_n-Function for n = 1
; http://mathworld.wolfram.com/En-Function.html
; It defined as special case of the exponential integral Ei(x)
; E1(x) = -Ei(-x)
; Rewritten from ExponentialIntegralEN (ALGLIB)
Func E1($x)
Local $Result, $r, $t
Local $yk, $xk, $pk, $pkm1, $pkm2, $qk, $qkm1, $qkm2
Local $psi = 0, $z, $i, $k, $big, $EUL
Local Const $MachineEpsilon = 1.11e-16

    $EUL = 0.57721566490153286060
    $big = 1.44115188075855872e17
    If $x<0 Or $x>170 Or $x=0 then return -1

    If $x<=1 then
        $psi = -$EUL-Log($x)
        $z = -$x
        $xk = 0
        $yk = 1
        $pk = 0
        $Result = 0.0
        Do
         $xk += 1
         $yk *= $z/$xk
         $pk += 1
         If $pk<>0 then $Result += $yk/$pk
         If $Result<>0 then
          $t = Abs($yk/$Result)
         Else
          $t = 1
         EndIf
        Until $t<$MachineEpsilon

        $t = 1
        $Result = $psi*$t-$Result
        return $Result
    Else
        $k = 1
        $pkm2 = 1
        $qkm2 = $x
        $pkm1 = 1.0
        $qkm1 = $x+1
        $Result = $pkm1/$qkm1
        Do
            $k += 1
            If Mod($k,2)=1 then
              $yk = 1
              $xk = 1+($k-1)/2
            Else
              $yk = $x
              $xk = $k/2
            Endif
            $pk = $pkm1*$yk+$pkm2*$xk
            $qk = $qkm1*$yk+$qkm2*$xk
            If $qk<>0 then
                $r = $pk/$qk
                $t = Abs(($Result-$r)/$r)
                $Result = $r
            Else
                $t = 1
            EndIf
            $pkm2 = $pkm1
            $pkm1 = $pk
            $qkm2 = $qkm1
            $qkm1 = $qk
            If Abs($pk) > $big then
                $pkm2 /= $big
                $pkm1 /= $big
                $qkm2 /= $big
                $qkm1 /= $big
            EndIf
        Until $t<$MachineEpsilon
        $Result *= exp(-$x)
    EndIf
    return $Result
EndFunc