Fixed Point implementation of C exponent function

[banu]

Hi,
I tried to google a lot searching for fixed point implementation of
exponent function y = exp(x) (e to the power x) but could not find it.
Also I think there is no post on this group also. I found some CORDIC
implementations but those were in C++. I need lightweight
implementation in C.

I want to use this function in an image processing algorithm to be
implemented on a VLIW-DSP which doesn't have a floating point unit.

I will appreciate in anyone here can suggest a plain C - code which
handles both positive and negative arguments (x<0 and x<=0) ,with
reasonable accuracy and performance.

thanks
varun

 

[Lew Pitcher]

Hi,
I tried to google a lot searching for fixed point implementation of
exponent function y = exp(x) (e to the power x) but could not find it. [snip]
I want to use this function in an image processing algorithm to be
implemented on a VLIW-DSP which doesn't have a floating point unit.

[snip]

When you say "fixed point", do you mean "short int / int / long int"?

ISTM that, because e is an irrational number (2.718281828459045...),
you will be in for a lot of inaccuracies if you are looking for an
integer implementation of exp().

OTOH, if an integer implementation is acceptable, then you are just
looking for a limited "positive power of two" sort of implementation
(ints can't go fractional, so -ve powers would result in zero).

 

[Ben Bacarisse]

I tried to google a lot searching for fixed point implementation of
exponent function y = exp(x) (e to the power x) but could not find it.
Also I think there is no post on this group also. I found some CORDIC
implementations but those were in C++. I need lightweight
implementation in C.

I think the answer will depend on the fixed-point format. If the
number of bits in the integer part is small (relative to available
memory), you can store the values of exp(i) for all of these and
re-write your exp(x) as exp(i)*exp(f). You can always choose i so
that f is < 0.5 which means that the normal Taylor series for exp(f)
will converge rapidly.

If you have a large range for x (relative to available memory) this
won't work and I am not sure what the best way to proceed is.

Either way, I have not given you what you want which is C code. Sorry
about that. I am sure there is some out there somewhere!

<snip>

 

[Ben Bacarisse]

I think the answer will depend on the fixed-point format. If the
number of bits in the integer part is small (relative to available
memory), you can store the values of exp(i) for all of these and
re-write your exp(x) as exp(i)*exp(f).

Correction: I think you will always be able to do this. Even with,
say, 64 bits in the integer part you will have to store fewer than 64
values of exp(i) because e > 2.

<snip>

 

[BGB / cr88192]

Hi,
I tried to google a lot searching for fixed point implementation of
exponent function y = exp(x) (e to the power x) but could not find it.
Also I think there is no post on this group also. I found some CORDIC
implementations but those were in C++. I need lightweight
implementation in C.

I want to use this function in an image processing algorithm to be
implemented on a VLIW-DSP which doesn't have a floating point unit.

I will appreciate in anyone here can suggest a plain C - code which
handles both positive and negative arguments (x<0 and x<=0) ,with
reasonable accuracy and performance.

well, as others have noted, it depends a lot on the number of bits in the
integer and fractional parts.

one possible idea that comes to mind is this:
rebase the exponent into 2^x form (this involves multiplying x by a constant
of 1/ln(2) or about 1.4426950...).

this way, a base can be calculated which is simply an integer power of 2.

then, one only need calculate a scale based on the fractional part, and if
the number of fractional bits is sufficiently small (or memory sufficiently
large), then a lookup table can be used.

then, the 2 numbers can be multiplied together, and one is done.

granted, for larger exponents this strategy may become less accurate. could
be addressed by using a bigger lookup table, or by using a recursive
strategy to calculate the fraction. similarly, some recursion could also
allow a smaller table at the cost of some speed.

another had mentioned using the taylor series, which is also an option, but
not likely as fast as using a lookup table.


for example (hypothetical):
#define FIX12_INVLN2 5909 //1.442695, 12-bit fraction
#define FIX28_INVLN2 387270501 //1.442695, 28-bit fraction

typedef int32_t fix12;
typedef int32_t fix28;

fix12 fix12_mul(fix12 a, fix12 b)
{
return ((fix12)((((int64_t)a)*b)>>12));
}

fix12 fix12_mul28(fix12 a, fix28 b)
{
return ((fix12)((((int64_t)a)*b)>>28));
}

static fix28 fix12_exp2lut[4096]={...}; //takes ~16kB memory

fix12 fix12_exp2(fix12 x)
{
fix12 b, f;
b=1<<(12+(x>>12)); //calculate base
f=fix12_exp2lut[x&4095];
return(fix12_mul28(b, f));
}

fix12 fix12_exp(fix12 x)
{ return(fix12_exp2(fix12_mul28(x, FIX28_INVLN2))); }


note:
some of the usual optimizations for optimizing fixed point calculations (or
improving accuracy) could possibly also be used here.

 

[Thad Smith]

The first thought is that if the C++ version is otherwise adequate, simply
convert the code to C.

well, as others have noted, it depends a lot on the number of bits in the
integer and fractional parts.
Agreed.

one possible idea that comes to mind is this:
rebase the exponent into 2^x form (this involves multiplying x by a constant
of 1/ln(2) or about 1.4426950...).

this way, a base can be calculated which is simply an integer power of 2.

That sounds good.

then, one only need calculate a scale based on the fractional part, and if
the number of fractional bits is sufficiently small (or memory sufficiently
large), then a lookup table can be used.

then, the 2 numbers can be multiplied together, and one is done.

granted, for larger exponents this strategy may become less accurate. could
be addressed by using a bigger lookup table, or by using a recursive
strategy to calculate the fraction. similarly, some recursion could also
allow a smaller table at the cost of some speed.

another had mentioned using the taylor series, which is also an option, but
not likely as fast as using a lookup table.


for example (hypothetical):
#define FIX12_INVLN2 5909 //1.442695, 12-bit fraction
#define FIX28_INVLN2 387270501 //1.442695, 28-bit fraction

typedef int32_t fix12;
typedef int32_t fix28;

fix12 fix12_mul(fix12 a, fix12 b)
{
return ((fix12)((((int64_t)a)*b)>>12));
}

fix12 fix12_mul28(fix12 a, fix28 b)
{
return ((fix12)((((int64_t)a)*b)>>28));
}

static fix28 fix12_exp2lut[4096]={...}; //takes ~16kB memory

fix12 fix12_exp2(fix12 x)
{
fix12 b, f;
b=1<<(12+(x>>12)); //calculate base
f=fix12_exp2lut[x&4095];
return(fix12_mul28(b, f));
}

fix12 fix12_exp(fix12 x)
{ return(fix12_exp2(fix12_mul28(x, FIX28_INVLN2))); }

I suggest the following changes to this reasonable approach:

The OP asked for a lightweight implementation, implying that he is attempting to
minimize code space and maybe data space.

The lookup table should be const. That avoids, in many implementations, having
separate RAM variable and initialization table. Since the values of minimum
differences between exp2lut entries is less than one in 6000, 28 fractional bits
is overkill. A 4.12 (or 1.15) representation in a short int cuts the table size
in half.

If code space is more important that execution time, two tables with for 6 bits
each (128 table entries) or three table entries of 4 bits each (48 total
entries) would reduce the constant data area. For example, have three tables:
2^(0.xxxx), 2^(0.0000xxxx), and 2^(0.00000000xxxx) where each position is a
single bit and ^ indicates exponentiation. Multiply the three relevant entries.

 

[BGB / cr88192]

The first thought is that if the C++ version is otherwise adequate,
simply convert the code to C.

yep, should work fine so long as there are not piles of classes or templates
in the mix...

That sounds good.

I suggest the following changes to this reasonable approach:

The OP asked for a lightweight implementation, implying that he is
attempting to minimize code space and maybe data space.

The lookup table should be const. That avoids, in many implementations,
having separate RAM variable and initialization table. Since the values
of minimum differences between exp2lut entries is less than one in 6000,
28 fractional bits is overkill. A 4.12 (or 1.15) representation in a
short int cuts the table size in half.

granted, this is an idea...

I guess 28 bits is overkill if there are only 4096 entries, but really, this
was intended as an example rather than an actually likely implementation
(after all, picking 12 fractional bits on my part was, really, rather
arbitrary...).

If code space is more important that execution time, two tables with for 6
bits each (128 table entries) or three table entries of 4 bits each (48
total entries) would reduce the constant data area. For example, have
three tables:
2^(0.xxxx), 2^(0.0000xxxx), and 2^(0.00000000xxxx) where each position is
a single bit and ^ indicates exponentiation. Multiply the three relevant
entries.

yep, also possible...

 

[Nick Keighley]

When you say "fixed point", do you mean "short int / int / long int"?

I guess he means "fixed point". This is a format for representing a
number with a certain number of bits before the decimal pojnt and a
certain number after. The decimal point is "fixed" as opposed to
"floating point" where it can move about. Fixed point arithmatic is
quite respectable and is often used on small embedded systems.

ISTM that, because e is an irrational number (2.718281828459045...),
you will be in for a lot of inaccuracies if you are looking for an
integer implementation of exp().

by that logic anything involving pi can't be calculated accuratly.
These things *can* be calculated to resonable accuracy in many
situations. Otherwise computer games wouldn't work.

OTOH, if an integer implementation is acceptable, then you are just
looking for a limited "positive power of two" sort of implementation
(ints can't go fractional, so -ve powers would result in zero).

this may not be what he wants.

 

[Malcolm McLean]

Hi,
I tried to google a lot searching for fixed point implementation of
exponent function y = exp(x) (e to the power x) but could not find it.

Go onto my website
http://www.personal.leeds.ac.uk/~bgy1mm

Go to Basic Algorithms, sample chapters, and look up the free chapter,
"maths functions". It tells you how to implement the standard library
maths functions. (It's not meant to be cut and paste code for
practical use, but example code to illustrate the general principles.
Sometimes extreme inputs are handled incorrectly, for instance. You
ough to be able to write your own fixed exp function using the exaple
as a guide).

 

[thadula]

I read your post when looking for fixed point exp(x) implementation. Eventually I ended up implementing it myself based on Texas Instruments guidelines. I know it's late response and most likely you already had your answer long ago, but just in case I'm posting my implementation of basic functions in 1.15 / 17.15 fixed point arithmetics. As some of those functions (that I required to be extremely fast) are based on lookup tables, using them requires initialization: fixedpoint_init(). The whole thing conatins implementation of sin/cos, atan, sqrt, ln, exp and some other simpler functions. As the post doesn't tolerate .h and .c attachments, I'm attaching it as .txt files.

Best regards,
Tom Hadula

 

[mcdivyakanth]

Hello Tom Hadula,

I would like to ask you if you can explain me/ share with me on how the look-up table for fraction part is arrived? ( for "fractionexptable" and for fractionexptable)
Thanks.