I am trying to write an algorithm to accurately calculate exponents (antilogs) for a variable precision floating point library I am working on. The base is not relevant since I can convert between them.

I was able to manually calculate log10() using repetitive application of x^10. This is a digit by digit calculation and requires 4 multiplies per digit. I can reverse the algorithm to calculate exp10(), but this requires repeated application of a 10th root. Calculating the 10th root is significantly more CPU costly than 10th power.

I searched the web and a lot of people suggested using a Taylor Series to calculate exp_e(). I did that and found that it requires about 2 iterations per digit for accurate results. Only two multiplies and one divide per iteration. This is still a bit steep in terms of CPU cycles especially when some FP numbers can be 100 digits long.

Now, I also found the algorithm that was used to calculate EXP in the old Sinclair ZX81. The author claimed that it was Chebyshev polynomials. I mention this because when I tested it, the algorithm was calculating accurately to one digit per iteration – much better than the Taylor Series.

I would use the algorithm as-is if it weren’t for the fact that the floating point library has to be accurate to an arbitrary number of digits. The ZX81 EXP code is only accurate to 8 digits. There is no explanation as to how to extend the number of iterations to get more accuracy.

So does anyone know how to calculate EXP() using Chebyshev Polynomials? Can they be expanded like the Taylor Series for more accuracy? Anything better than either?

(Please no long math proofs. That’s over my head. I just want the algorithms.)