I have extended the algorithm for the Exponential Integrals with a routine to
evaluate the functions for Bigfloat numbers. The algorithm seems to work very
fine for real and complex Bigfloat numbers. But I have problems to verify the
results.
I have taken examples from functions.wolfram.com with a precision of 64 digits.
But I can only reproduce the values in the best case within an accuracy of 31
digits and in the worst case within an accurarcy of 16 digits.
Here an example for the best case with E1(1.5):
At first Maxima double float calculation:
0.10001958240663
Next Maxima Bigfloat with 16,32 and 64 digits (continued fractions):
1.000195824066326b
1.0001958240663265190190933991165b-1
1.000195824066326519019093399116400594658618543362758786084650014b-1
functions.wolfram.com with 64 digits:
0.1000195824066326519019093399116669782617300061403505850505670604
The results of Maxima agree within 31 digits.
Here an example for the worst case with E1(0.5):
At first Maxima double float calculation:
0.55977359477616
Next Maxima Bigfloat with 16,32 and 64 digits (series expansion):
5.597735947761609b-1
5.5977359477616080680388078688444b-1
5.597735947761608068038807868844401350218672492054958640540604543b-1
functions.wolfram.com with 64 digits:
0.5597735947761608117467959393150852352268468903163535152482932191
The results of Maxima agree within 16 digits.
I don't know the algorithm of functions.wolfram.com, but I hope they give
numerical values which are correct.
I suppose the problems arise because we need the constants BIGFLOAT%E and
BIGFLOAT%GAMMA in the calculation. By inspection of the file float.lisp, I have
seen that these constants are implemented with an accuracy of 16 digits. The
limited accuracy of these constants might be the bottle neck of the algorithm.
I have no experience to calculate the mantissa of Bigfloat numbers. Perhaps,
someone could help me to get these constants with much higher accuracy as
Bigfloats.
Dieter Kaiser