\ fexpint     Real Exponential Integral         ACM Algorithm #20

\ Forth Scientific Library Algorithm #1

\ Evaluates the Real Exponential Integral,
\     E1(x) = - Ei(-x) =   int_x^\infty exp^{-u}/u du      for x > 0
\ using a rational approximation

\ This code conforms with ANS requiring:
\      1. The Floating-Point word set

\ Collected Algorithms from ACM, Volume 1 Algorithms 1-220,
\ 1980; Association for Computing Machinery Inc., New York,
\ ISBN 0-89791-017-6

\ (c) Copyright 1994 Everett F. Carter.  Permission is granted by the
\ author to use this software for any application provided the
\ copyright notice is preserved.

\ CR .( EXPINT     V1.1                  21 September 1994   EFC )
\ include lib/ansfloat.4th

[UNDEFINED] fexpint [IF]
[UNDEFINED] fexp [IF] include lib/fexp.4th    [THEN]
[UNDEFINED] fln  [IF] include lib/flnflog.4th [THEN]
: FEXPINT ( --, f: x -- expint[x] )
        FDUP
        1 s>f F< IF
                    FDUP s" 0.00107857" s>float F* s" 0.00976004" s>float F-
                    FOVER F*
                    s" 0.05519968" s>float F+
                    FOVER F*
                    s" 0.24991055" s>float F-
                    FOVER F*
                    s" 0.99999193" s>float F+
                    FOVER F*
                    s" 0.57721566" s>float F-
                    FSWAP FLN F-
                ELSE
                    FDUP s" 8.5733287401" s>float F+
                    FOVER F*
                    s" 18.059016973" s>float F+
                    FOVER F*
                    s" 8.6347608925" s>float F+
                    FOVER F*
                    s" 0.2677737343" s>float F+

                    FOVER
                    FDUP s" 9.5733223454" s>float F+
                    FOVER F*
                    s" 25.6329561486" s>float F+
                    FOVER F*
                    s" 21.0996530827" s>float F+
                    FOVER F*
                    s" 3.9584969228" s>float F+

                    FSWAP FDROP
                    F/
                    FOVER F/
                    FSWAP -1 s>f F* FEXP
                    F*

                THEN
;

\ test code generates a small table of selected E1 values.
\ most comparison values are from Abramowitz & Stegun,
\ Handbook of Mathematical Functions, Table 5.1

false [IF]
: fexpint_test ( -- )
      CR
      ."  x     E1(x) exact    x FEXPINT " CR

      ."  0.5   0.5597736      "
      s" 0.5" s>float fexpint  F. CR

      ."  1.0   0.2193839      "
      1 s>f fexpint   F. CR

      ."  2.0   0.0489005      "
      2 s>f fexpint    F. CR

      ."  5.0   0.001148296    "
      5 s>f  fexpint    F. CR

      ." 10.0   0.4156969e-5   "
      10 s>f fexpint   FS. CR
;

fclear
10 set-precision
fexpint_test
[THEN]
[THEN]