Fixed Point implementation of C exponent function

Discussion in 'C Programming' started by banu, Mar 7, 2010.

  1. banu

    banu Guest

    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
     
    banu, Mar 7, 2010
    #1
    1. Advertising

  2. banu

    Lew Pitcher Guest

    On Mar 7, 9:33 am, banu <> wrote:
    > 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).
     
    Lew Pitcher, Mar 7, 2010
    #2
    1. Advertising

  3. banu <> writes:

    > 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.
     
    Ben Bacarisse, Mar 7, 2010
    #3
  4. Ben Bacarisse <> writes:

    > banu <> writes:
    >
    >> 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).


    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>
    --
    Ben.
     
    Ben Bacarisse, Mar 7, 2010
    #4
  5. "banu" <> wrote in message
    news:...
    > 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.


    > thanks
    > varun
     
    BGB / cr88192, Mar 7, 2010
    #5
  6. banu

    Thad Smith Guest

    BGB / cr88192 wrote:
    > "banu" <> wrote in message
    > news:...
    >> 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.


    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.

    --
    Thad
     
    Thad Smith, Mar 8, 2010
    #6
  7. "Thad Smith" <> wrote in message
    news:4b945e05$0$65843$...
    > BGB / cr88192 wrote:
    >> "banu" <> wrote in message
    >> news:...
    >>> 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.

    >
    > 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...


    >> 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.
    >


    <snip>

    >
    > 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...


    > --
    > Thad
     
    BGB / cr88192, Mar 8, 2010
    #7
  8. On 7 Mar, 15:31, Lew Pitcher <> wrote:
    > On Mar 7, 9:33 am, banu <> wrote:>


    > > 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..

    >
    > > 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.

    >
    > 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.
     
    Nick Keighley, Mar 8, 2010
    #8
  9. On 7 Mar, 14:33, banu <> wrote:
    > 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).
     
    Malcolm McLean, Mar 8, 2010
    #9
  10. banu

    thadula

    Joined:
    Dec 24, 2010
    Messages:
    1
    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
     

    Attached Files:

    thadula, Dec 24, 2010
    #10
  11. banu

    mcdivyakanth

    Joined:
    Oct 15, 2012
    Messages:
    1
    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.
     
    mcdivyakanth, Nov 29, 2012
    #11
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. H aka N
    Replies:
    15
    Views:
    15,724
    Ben Jones
    Mar 2, 2006
  2. Motaz Saad
    Replies:
    7
    Views:
    6,516
  3. Replies:
    6
    Views:
    401
  4. Saraswati lakki
    Replies:
    0
    Views:
    1,374
    Saraswati lakki
    Jan 6, 2012
  5. Marco
    Replies:
    2
    Views:
    111
    Serhiy Storchaka
    Jul 16, 2013
Loading...

Share This Page