/**
   Program to calculate the log integral, often called Li[x]

   This is Integrate[1/ln[x], 0, x]

   This goes to -infinity at 1 so don't be surprised if it explodes there.

   This is useful in number theory, such as estimating the number of primes
   below a certain number.

   This function does not currently use arbitrary precision for ln[z] so
   its precision is limited and also limited by the precision of Frink's
   internally-stored value of EulerMascheroniConstant.  But this version is
   relatively fast for about 14 digits of precision.

   If you need arbitrary precision, it will be better to use Frink's
   arbitrary-precision numerical integration, for example:

   use NumericalIntegrationArbitrary.frink
   digits = 80
   setPrecision[digits]
   LogIntegral[a,b,digits=getPrecision[] ] := NIntegrate[{|z| 1/arbitraryLn[z]}, a, b, 10^-digits]
   LogIntegral[2, 10, digits]

   Currently, it will be necessary to set "digits" to a few more than the
   desired number of digits to make the last few digits come out right.

   See:  http://functions.wolfram.com/GammaBetaErf/LogIntegral/
*/
LogIntegral[z] :=
{
   lnz = ln[z]
   sum = (1/2) (ln[lnz] - ln[1/lnz]) + EulerMascheroniConstant

   k = 0
   
   do
   {
      k = k+1
      oldsum = sum
      term = lnz^k / (k k!)
      sum = sum + term
      //  println["$k\t$term\t$sum"]
   } while (oldsum != sum)

   return sum
}