Computing $e^x$ in Octave using only +,-,*,/,^.

Question

I have to write a function in Octave which is computing e^x. I wrote this:

function y=f(x)
    y=0; t=1; k=1;
        while abs(t)>0.001
        y=y+t;
        t=t*(x/(2*k));
        k=k+1;
        end
    y=y^2;
end;

Everything is okay when I'm computing $f(20)$ or $f(-30)$ etc. But when I try to compute $f(-2000)$ I should get $ \approx 0$, but Octave says it's $\infty$. Is it my fault (I mean something's wrong with code)?

The problem is floating point numbers are not real numbers. And for large negative numbers, there is much cancelation in this code. See http://www.validlab.com/goldberg/paper.pdf — Jean-Claude Arbaut, Mar 23 '13 at 22:21
have you tried debugging? (like printing the intermediates each loop) — V-X, Mar 23 '13 at 22:23
To compute $e^x$ with taylor expansion, manage to use only small $x$. To achieve this, use $e^x=\left(e^{x\over n}\right)^n$, so that argument in power series may be kept small. You allready use this $n=2$, but you have to chose a $n$ near $x$. — Jean-Claude Arbaut, Mar 23 '13 at 22:23
Also consider computing $e^{-x} = 1 / e^x$ to avoid cancellations. Or use a CORDIC algorithm for kicks... — vonbrand, Mar 23 '13 at 22:47

score 1 · Answer 1 · answered Mar 24 '13 at 02:34

In the spirit of arbautjc's comment, given the operators you have, you can always ensure that the magnitude of the argument is less than $1$. Furthermore if $x<0$, you can do as vonbrand suggested and use index laws to aid your calculations. This will likely be much slower than simply using the exp(..) command built in.

The following is MATLAB code, but it should work in Octave as well.

function y = myexp(x)
    % determine if x is negative or not
    if (x < 0)
        isneg = true;
        x = -x;
    else
        isneg = false;
    end

    % scale the input variable to have magnitude less than 1
    scale = 1; % initialise the scale argument
    factor = 16; % set a factor which will be used to increment the scale
    while x > 1
        x = x / factor;
        scale = scale * factor;
    end

    % evaluate exp(x) with 0<x<1
    y=0; t=1; k=1; % initialise parameters
    tol = 1e-6; % convergence tolerance
    while abs(t)>tol
        y=y+t;
        t=t*(x/k);
        k=k+1;
    end

    % unscale the result
    y=y^scale;

    % check if the input was negative and reciprocate the result if needed
    if isneg
        y = 1/y;
    end
end

Some results are below. V is value and T is execution time. tol was set to 1e-6.

x=-100
exp V: 3.720076e-44 T: 0.000004
myexp   V: 3.720262e-44 T: 0.000059
x=-10
exp V: 4.539993e-05 T: 0.000002
myexp   V: 4.540017e-05 T: 0.000021
x=-5
exp V: 6.737947e-03 T: 0.000002
myexp   V: 6.737952e-03 T: 0.000016
x=-4
exp V: 1.831564e-02 T: 0.000002
myexp   V: 1.831572e-02 T: 0.000016
x=-3
exp V: 4.978707e-02 T: 0.000002
myexp   V: 4.978711e-02 T: 0.000016
x=-2
exp V: 1.353353e-01 T: 0.000002
myexp   V: 1.353358e-01 T: 0.000016
x=-1
exp V: 3.678794e-01 T: 0.000002
myexp   V: 3.678795e-01 T: 0.000015
x=0
exp V: 1.000000e+00 T: 0.000002
myexp   V: 1.000000e+00 T: 0.000016
x=1
exp V: 2.718282e+00 T: 0.000002
myexp   V: 2.718282e+00 T: 0.000016
x=2
exp V: 7.389056e+00 T: 0.000002
myexp   V: 7.389029e+00 T: 0.000016
x=3
exp V: 2.008554e+01 T: 0.000002
myexp   V: 2.008552e+01 T: 0.000016
x=4
exp V: 5.459815e+01 T: 0.000002
myexp   V: 5.459791e+01 T: 0.000017
x=5
exp V: 1.484132e+02 T: 0.000002
myexp   V: 1.484131e+02 T: 0.000016
x=10
exp V: 2.202647e+04 T: 0.000002
myexp   V: 2.202635e+04 T: 0.000017
x=100
exp V: 2.688117e+43 T: 0.000002
myexp   V: 2.687982e+43 T: 0.000016

Computing $e^x$ in Octave using only +,-,*,/,^.

1 Answers1