Update FermiDiracOneHalfTeSCA

Replace evaluations of log(1+x) with log1p(x).
In the code this occurs twice.

In the old version, trickery in the first branch catched the case
when 1+x == 1 with a test if w==0.

We had this discussion before, the trick came from
JuliaDiff/ForwardDiff.jl#481

This does not work anymore with
JuliaDiff/ForwardDiff.jl#481

This PR (unfortunately) is now in v0.10.33 of ForwardDiff.jl.
There is a discussion going on that this should have been 0.11, i.e. marked
as a breaking change.

It did break Example 103, the PR fixes this.

Author: j-fu

src/ct_distributions.jl

@@ -1,4 +1,3 @@
-
 using ForwardDiff
 
 """
@@ -65,24 +64,23 @@ $(TYPEDSIGNATURES)
 
 The incomplete Fermi-Dirac integral of order 1/2, implemented according to the software
 package TeSCA, see https://wias-berlin.de/software/index.jsp?lang=1&id=TeSCA.
+
+Modified to use log1p(x)=log(1+x).
 """
 function FermiDiracOneHalfTeSCA(x::Real)
     if x < 1.6107
-        ex = exp(x)
-        y  = 1+ex
-        w  = y-1
-        z  = w==0 ? ex : ex*log(y)/w
+        z=log1p(exp(x))
         return ( 1 + 0.16 * z ) * z
     elseif 1.6107 <= x <= 344.7
-        z = log( 1 + exp( x^(3/4)) )
-        return 0.3258 - 0.0321 * z  + 0.7523 * z^2
+        z = log1p(exp( x^(3/4)) )
+        return 0.3258 - (0.0321  - 0.7523 * z ) * z
     else
         z = x^(3/4)
-        return 0.3258 - 0.0321 * z  + 0.7523 * z^2
+        return 0.3258 - (0.0321 -  0.7523 * z ) * z
     end
-
 end
 
+
 """
 $(TYPEDSIGNATURES)