How to compute exp(x) for x in [0,ln(2)] in double precision with high accuracy, possibly at most 0.5 ulp error?

taumulurtulkyoy 2022-10-23 Answered
How to compute exp(x) with correct rounding?
How to compute exp(x) for x∈[0,ln(2)] in double precision with high accuracy, possibly at most 0.5 ulp error?
It is possible to obtain low error. After testing a few libraries I have found one with 1 ulp accuracy. However I am not able to find any reference, which delivers precise error bounds for computed value.
Calculating exp(x) is not a trivial task. Althougth is is easy to obtain fast converging polynomial approximating this function, for example the Taylor series, roundoff errors occuring during evaluation of a polynomial are quite high. For example if we want to evaluate a polynomial p(x)=∑N−1i=0aixi, then apriori relative forward error is bounded as
| fl ( p ( x ) ) p ( x ) p ( x ) | cond ( p , x ) × u + O ( u 2 )
where
cond ( p , x ) = 2 N 1 2 N u × | i = 0 N | a i | | x | i p ( x ) | 2 N
and u is the unit roundoff (half of the machine precision).
Therefore even if we have a polynomial approximation of exp(x) of order 5 (which is hard by itself), with truncation error at most 0.5 ulp, then still roundoff error is at least 5. However 1 ulp error is possible.
My question is what kind of tricks can be used here to obtain high accuracy (without using high precision arithmetics). Is it even possible to obtain 0.5 ulp error?
You can still ask an expert for help

Want to know more about Polynomial arithmetic?

Expert Community at Your Service

  • Live experts 24/7
  • Questions are typically answered in as fast as 30 minutes
  • Personalized clear answers
Learn more

Solve your problem for the price of one coffee

  • Available 24/7
  • Math expert for every subject
  • Pay only if we can solve it
Ask Question

Answers (1)

Messiah Trevino
Answered 2022-10-24 Author has 18 answers
Step 1
I have found one nice formula for the exp function:
(1) G ( x ) = coth ( x / 2 ) x 2 x (2) c = x G ( x 2 ) × x 2 (3) exp ( x ) = 1 + ( x c 2 c + x )
This formula was used for interval [-ln(2)/2,ln(2)/2] but this is also good (transformation from this interval to [0,ln(2)] can be done without any errors. Unfortunately I don't have any reference on this scheme. I found this scheme using "reverse engeneering" of a library source.
The function G(x) is nearly linear in this interval. This function can be accurately approximated by a polynomial of order 4 using the Remez algorithm. Its Taylor expansion is
G ( x ) 1 6 x 360 + x 2 15120 x 3 604800 + x 4 23950080
Step 2
Condition number of this polynomial in given interval is very good, 1.012 × γ 5 , γ N = 2 N / ( 1 2 N u ).
Expression (2) cannot be computed accurately. However (3) suppresses most of numerical errors. For example condition number of (3) with respect to c is 0.0619, thus error if computing c, δ c , contributes to the total error as 0.0619 × δ c , similarly x c / ( 2 c ) is inaccurate, but error in computing this subexpression δ x c contributes to error of x c / ( 2 c ) + x approximately as 0.2 × δ x c , which again is suppressed in the final sum (3) (approximately twice).
Final error of computing exp by (1)-(3) is less than 1 (maximum error found is 0.83).
It is possible to compute exp(x) even with 0.5 ulp error using the compensated Horner's scheme. See "Faithful Polynomial Evaluation with Compensated Horner Algorithm", P. Langlois, N. Louvet, 2006, for details. This scheme must however be corrected in order to take into account errors in computed coefficients of approximation polynomial. As the polynomial approximation of exp(x) we can use the Taylor expansion of order 27. Probably we can use lower order approximation, in tests I was not able to find a case with higher than 0.5 ulp even in case of approximation of order 20. Moreover this scheme is quite fast, at most 8-times slower, than library implementation of exp.
Did you like this example?
Subscribe for all access

Expert Community at Your Service

  • Live experts 24/7
  • Questions are typically answered in as fast as 30 minutes
  • Personalized clear answers
Learn more

You might be interested in

asked 2021-03-02
Write A if the sequence is arithmetic, G if it is geometric, H if it is harmonic, F if Fibonacci, and O if it is not one of the mentioned types. Show your Solution. a. 13,29,327,481,... b. 3,8,13,18,...,48
asked 2021-08-14

Find the total number of terms in the given arithmetic sequence.
a1=296,d=13,S=36

asked 2020-11-23
1. Explain with numerical examples what Real Numbers and Algebraic Expressions are. 2. Explain with numerical examples Factoring and finding LCMs (least common multiples). Explain factoring of a larger number. 3. Explain with numerical examples arithmetical operations (addition, subtraction, multiplication, division) with fractions 4, Explain with numerical examples arithmetical operations (addition, subtraction, multiplication, division) with percentages 5. Explain with numerical examples exponential notation 6. Explain with numerical examples order (precedence) of arithmetic operations 7. Explain with numerical examples the concept and how to find perimeter, area, volume, and circumference (use related formulas)
asked 2021-02-20

Indicate whether the expression defines a polynomial function. P(x)=x2+3x+3 polynomial or not a polynomial If it is a polynomial function, identify the following. (If it is not a polynomial function, enter DNE for all three answers.)

(a) Identify the leading coefficient.

(b) Identify the constant term.

(c) State the degree.

asked 2022-05-08
Solve 21-?=19.
asked 2021-07-29
To calculate: The coefficient and the degree of each term of the polynomial x2y+4y32xy and then find the degree of the polynomial.
asked 2022-11-06
Some search on the internet and this site didn't result in any topic about this question of Silverman's The Arithmetic of Elliptic Curves:
Let W P n be a smooth algebraic set, each of whose component varieties has dimension n 1. Prove that W is a variety.
Hint: use Krull's Hauptidealsatz to show W is the zeroset of a single homogenous polynomial.
my ideas: So if I(V) contains at least two linear independent polynomials, we have to prove the height of I(V) is at least 2.
After that, it is enough to find one singular point, where we assume I ( V ) = f g with f,g nonconstant polynomials.