Finding the exponential of a very small number

Hi

I wish to find the exact exponential of numbers in the range of 10^-23. Matlab gives me the answer as 1 for exp(3.528393650226818e-23). This is not right. How do I accomplish this?

Answer 1

[John D'Errico]

Multiple ways to solve this problem. Sorry, but quad precision is not an option. There is no way to achieve quad precision, at least without recourse to some other tool, or exporting the results to a language where quad precision is available.

Brute force, without resource to the symbolic toolbox ... Use my HPF toolbox, which is free to download. I'll use 100 decimal digits of precision, even though

DefaultNumberOfDigits 100
s = hpf('-6.471606349773182e-23');
d = hpf('-6.604939683106516e-23');

See that HPF has no problem in evaluating expressions like those that you have.

exp(s)
ans =
0.9999999999999999999999352839365022681800000020940844373212284440930971466223949024568510358194049670

log(1+exp(s))
ans =
0.6931471805599453094171997634264277021655001348837763634509871205166792020337447156058633269050605702

So A is just

A = hpf('1.449056603773585e+19')*(log(1+exp(s))-log(1+exp(d)))
A =
0.000009660377358490614968553143308576941771731705532455911053800738735222612516743226005743174365407550000

As a check with the symbolic toolbox, it yields:

s = sym('-6.471606349773182e-23');
d = sym('-6.604939683106516e-23');
A = sym('1.449056603773585e+19')*(log(1+exp(s))-log(1+exp(d)))
A =
0.0000096603773584906138365242807537049

In fact, the symbolic toolbox loses some precision here, because it used limited precision to do some of those computations. And, since you only provided those numbers to about 20 significant digits, what can you expect? But it does reasonably well. The exact result, reporting it as 32 significant digits is:

0.0000096603773584906149685531433085769

Can we get that result by using only numerical methods in double precision? Brute force will of course fail, as you found.

exp(-6.471606349773182e-23)
ans =
     1

A = 1.449056603773585e+19*(log(1+exp(s))-log(1+exp(d)))
A =
     0

So in order to get a viable result for A using double precision, one needs to resort to numerical analysis. We need to understand that even if we could compute this for very tiny numbers s and d, the results log(1+exp(s)) and log(1+exp(d)) are very close to each other. Using HPF, the exact results out to 40 decimal digits are:

log(1+exp(s))
ans =
0.6931471805599453094171997634264277021655

log(1+exp(d))
ans =
0.6931471805599453094171990967597610354955

So, when we subtract the two results, we will see massive subtractive cancellation. HPF is able to give the exact results here because I did the computations in sufficiently high precision, but then only reported the first 40 digits.

Going back to double precision, even tools in MATLAB like log1p must fail, because exp(s) returns exactly 1. As far as MATLAB is concerned, with doubles for s and d, these two results are identical.

log(1+exp(s)) == log(1+exp(d))
ans =
  logical
   1

Ok, so can we do something? Yes. With some mental effort. Note that is s is TINY, then exp(s) is essentially 1. And then we lose accuracy because of the massive subtractive cancellation, because both logs are returning essentially log(2).

What if we computed the result of

log(1+exp(x)) - log(2)

Of course, we need to do this computation accurately, for very tiny x. First see that

log(1+exp(x)) - log(2) = log((1+exp(x))/2)

When we do the subtraction, those -log(2) terms will exactly cancel each other out. But they enable us to do the computation in double precision. We will need to use a Taylor series. I'll use the symbolic toolbox here, because I'm too lazy to do the algebra on paper.

syms x
taylor(log((1+exp(x))/2))
ans =
-x^4/192 + x^2/8 + x/2

For tiny values of x of magnitude smaller than eps, we can truncate that series to one term. But we can see what happens if we use 2 terms anyway.

logapprox1 = @(x) x/2;
logapprox2 = @(x) x/2 + x.^2/8;

s = -6.471606349773182e-23;
d = -6.604939683106516e-23;

format long g

A = 1.449056603773585e+19*(logapprox1(s)-logapprox1(d))
A =
      9.66037735849058e-06

A = 1.449056603773585e+19*(logapprox2(s)-logapprox2(d))
A =
      9.66037735849058e-06

As you can see, the two results are both reasonably accurate. Remember, the exact result (using scientific notation) was:

      9.6603773584906149685531433085769e-6

So our double precision computation was correct out to almost the full precision available using doubles.

The trick was to use finesse. Brute force only works if you are willing and able to do the computation using the proper tools that will support brute force. As I showed, even the symbolic toolbox missed a few digits out there at the end.

[Walter Roberson]

sym('-6.471606349773182e-23') turns out to represent the values as some kind of internal binary fraction that is not as simple as moving the decimal place and using an appropriate power of 10 as denominator. I did some analysis of the representation error in a post a few months ago, but I would have to dig it out.

Answer 2

[Christoph F.]

> This is not right.

Is is as close to the correct answer as can be represented with the 64-bit floating point format that Matlab uses.

exp(3.528393650226818e-23) is approximately 1 + 3.528393650226818e-23 (Taylor series, terminated after the first term), which will also give the result "1".

If 64-bit floating point is not enough precision for your purpose, you should look for something that supports 80-bit or quadruple precision (128-bit) floating point data types. Or redefine your task (i.e. how many digits are required?), since a computer will never give you an exact value for a number that has no finite representation.

[Himanshi Rani]

Are there any methods to find quadruple precision in matlab?

[Christoph F.]

I do not know of any official support for quadruple precision floating point in Matlab.

However, Jan Simon suggested using variable precision arithmetics, but this requires the symbolic math toolbox (which I do not have available, so I cannot comment on how useful it is for solving your task).

Can you define how many digits of precision you require? Approximating exp(x) as 1+x for very small x is already quite precise, as the next term of the Taylor series, x^2/2, is really tiny if x is small to begin with (if x is ~1e-23, then x^2 is ~1e-46).

[Himanshi Rani]

I have to actually find the value of A to enough precision:

A= (1.449056603773585e+19)*(log(1+exp(s))-log(1+exp(d)))

where s= -6.471606349773182e-23 and d= -6.604939683106516e-23

By default i.e without using the symbolic toolbox or the method described by Jan Simon, I get A as 0

Answer 3

[Walter Roberson]

taylor log(1+exp(x)) to order 5 (so you get the ^4 terms). evaluate at symbol x1 and at symbol x2 and subtract the two -- which will cancel out a log(2) . Multiply by 1.449056603773585e+19 . factor to get

   -75471698113207552 * ( x1 - x2) * ( x1^3 + x1^2*x2 + x1*x2^2 - 24*x1 + x2^3 - 24*x2 - 96)

You could work with that, or you could note that in your range 24*x1 or 24*x2 are going to be noise compared to 96, so taking into account floating point accuracy you could reduce this to

   -75471698113207552 * ( x1 - x2) * ( - 96)

Multiply out the two constants...

7245283018867924992 * (x1 - x2)

For your s and d, this formula is off in at most the 16th decimal place.

The 5th order taylor series differs from a high precision answer starting at 2^(-70)

[Himanshi Rani]

I could not understand your answer. Is x1= exp(s) ans x2=exp(d)? Also, I do not understand how you get log(2)

[Walter Roberson]

x1 is s and x2 is d. 7245283018867924992 * (x1 - x2) is the formula for the complete expression for A

... I might have had a bracket in the wrong place. I will look more carefully after I get some sleep.

[Walter Roberson]

No, I was right. You have

A= (1.449056603773585e+19)*(log(1+exp(s))-log(1+exp(d)))

which is a constant times log(1+exp(s)) - log(1+exp(d))

As s approaches 0, exp(s) approaches 1, so you approach

log(1+1) - log(1+exp(d))

which gets you the log(2) . The same argument for exp(d) with small d, so you headed towards something that very closely approximates log(2) minus something else that very closely approximates log 2. When you taylor the general form log(1+exp(x)) near 0 to get a bit more accuracy, one of the elements in the sum is going to be that expression evaluated at exactly 0, giving log(2); that part of the sum is independent of s or d, so when you subtract the two taylor forms, the log(2) cancel out, leaving only the terms that follow.

taylor(log(1+exp(x)),x,'Order',3)

is x^2/8 + x/2 + log(2) . Do not use direct precision arguments with respect to the log(2) because we are going to cancel those out as described above, so we are left considering x^2/8 + x/2 as the effective value. When x is in the order of 1e-23, then we know that x^2/8 will be on the order of 1e-47 which is negligible compared to 1e-23/2 . We can therefore see that an order 2 approximation is enough: near 1e-23, log(1+exp(x)) is approximately x/2 + log(2) . Substitute in s and d, subtract, and we see that log(1+exp(s)) - log(1+exp(d)) is going to be approximately s/2 - d/2 . All of which is to be multiplied by (1.449056603773585e+19)

Answer 4

[Jan]

The more exact numerical output is - as expected:

1.000000000000000000000035283936502268180000000622478087548...

If you do not want to use VPA and symbolical solutions, you can apply some mathematics also. Look at the Taylor expansion:

exp(x) = 1 + x + x^2/2 + x^3/6 + ...

For x very near to 0.0 this is 1 + x. As Cristoph has explained already 1 + 5.52e-23 cannot be represented "accurately" with the 64 bit IEEE 754 doubles. But perhaps you can proceed with your calculations using the 2 variables:

c1 = 1
c2 = x

and consider, that you need the sum of both.

[Himanshi Rani]

I have used the taylors expansion and it is good. But I do not know how to use the vpa and symbolical solutions. Could you please explain because I have to do a lot of computations with these numbers other than findidng exponential?

[Jan]

Do you have the Symbolic Toolbox? Unfortunately I do not have it. For a suggestion it would be useful if you explain exactly, what you want to achieve.