#  --------------------------------------------------------------
#
   FDm12_CT := proc(x::numeric)
#
#  --------------------------------------------------------------
#
#  PURPOSE
#
#    This Maple 6 procedure Computes the Fermi-Dirac integral
#
#       Fj(x) = 1/GAMMMA(j+1)*int(t^j/(exp(t-x)+1),t=0..infinity)
#
#    of order j = -1/2, correct to at least 7 significant decimal
#    digits of accuracy using Cody and Thacher's piecewise 
#    rational approximation scheme. (See reference below.)
#
#    The magnitude of the relative error is at most 0.25x10^(-7)
#    on the x domain (-infinity,infinity). 
#
#  USAGE
#
#    y = FDm12_CT(x)
#
#  FILE NAME 
#
#    FDm12_CT.txt
#
#  DATE
#
#    October 3, 2000
#
#  Programmer
#
#    Dr. F.G. Lether
#    Dept. Math. 
#    Univ. Georgia
#    Athens, GA. 30602 USA
#
#    email : fglether@math.uga.edu
#
#  REFERENCE
#
#    Cody,W.J. and Thacher,H.C.(1967),Rational Chebyshev
#    Approximations for Fermi-Dirac Integrals of Orders -1/2,
#    1/2 and 3/2, Math. Comp. 21, 30-40.
#
#  --------------------------------------------------------------
#
   local const,p1,q1,p2,q2,p3,q3,R,t,FDm12;
   Digits := 13;
#
#  Initialize constant 1/GAMMA(1/2) = 1/sqrt(Pi)
#	
   const := 0.5641895835478;
#   
#  Region 1 : Left tail (-infinity,1] Cody-Thacher coefficients
#
   p1 := array(0..4 , [ -1.253314128820    ,
                        -1.723663557701    ,
                        -6.559045729258e-1 ,
                        -6.342283197682e-2 ,
                        -1.488383106116e-5 ] );
   q1 := array(0..4 , [  1.000000000000    ,
                         2.191780925980    ,
                         1.605812955406    ,
                         4.443669527481e-1 ,
                         3.624232288112e-2 ] );
#
#  Region 2 : Middle region [1,4] Cody-Thacher coefficients
#
   p2 := array(0..4 , [ 1.0738127694 ,
                        5.6003303660 ,
                        3.6882211270 ,
                        1.1743392816 ,
                        2.3641935527e-1 ] );
   q2 := array(0..4 , [ 1.0000000000    ,
                        4.6031840667    ,
                        4.3075910674e-1 ,
                        4.2151132145e-1 ,
                        1.1832601601e-2 ] );
#
#  Region 3 : Right tail [4,infinity) Cody-Thacher coefficients
#
   p3 := array(0..4 , [ -8.222559330e-1 ,
                        -3.620369345e1  ,
                        -3.015385410e3  ,
                        -7.049871579e4  ,
                        -5.698145924e4  ] );
   q3 := array(0..4 , [  1.000000000   ,
                         3.935689841e1 ,
                         3.568756266e3 ,
                         4.181893625e4 ,
                         3.385138907e5 ] );
#
#  Explicit function to evaluate Horner nested form of rational
#  function Rn,n(z) = Pn(x)/Qn(z), with n = 4
#
   R := (P,Q,z) ->
        ( P[0] + (P[1] + (P[2] + (P[3] + P[4]*z)*z)*z)*z )
      / ( Q[0] + (Q[1] + (Q[2] + (Q[3] + Q[4]*z)*z)*z)*z );
#
#  Compute Cody-Thacher approximation depending on x in
#  (-infinity,1), [1,4) or [4,infinity). (Use the Maple
#  hardware floating point evaluation instruction, evalhf,   
#  to gain speed.)
#
   if x < 1 then
     t := evalhf(exp(x));
     FDm12 := evalhf( t * ( 1 + const * t * R(p1,q1,t) ) );
   elif x < 4 then
     FDm12 := evalhf( const * R(p2,q2,x) );
   else
     t := evalhf( 1/(x^2) );
     FDm12 := evalhf( const * sqrt(x) * ( 2 + t * R(p3,q3,t) ) );
   fi;
#
#  Round floating point result to Digits = 13  significant digits,
#  with at least 7 of these significant digits being correct  
#
   FDm12 := evalf(FDm12,Digits);
#
   end: