johns@xetron.com (John Sampson) wrote in message news:<ae739955.0407070338.5482ba02@posting.google.com>...
> kenhusier@yahoo.com (Ken) wrote in message news:<46c498f7.0407060822.6197b9b7@posting.google.com>...
> > Mark Borgerding <mark@borgerding.net> wrote in message news:<MmmFc.6272$m91.5379@fe1.columbus.rr.com>...
> > > Ken wrote:
> > > > 
> > > > I have a question, if I need to calculate float point exp()function,
> > > I assume you mean fixed point exp(), right?
> > > 
> > > > what should I do? have to change them back to float, right? if then, I
> > > > still can't avoid from a floating point arithmatic.
> > > > 
> > > > Thanx,
> > > > Ken
> > > 
> > > You left out a VERY important piece of info.  What is your fixed point 
> > > representation?  I'm assuming Q15 , since implementing an exp() function 
> > > for fixed-point representation with more than one or two integer bits 
> > > becomes ridiculous.   exp(127) =~ 1e55
> > > 
> > > 
> > > First, decide how you'll represent the return value.
> > > 
> > > If you're using Q15 ints, the range is [-1:1).
> > > exp(-1) = 1/e =~ .367
> > > exp(1) = e = ~ 2.718
> > > 
> > > Assuming you can overcome a constant multiplier, I suggest returning 
> > > exp(x)/e to prevent overflow.
> > > 
> > > The output range becomes [ exp(-2):1 ] =~ [ .13534... : 1 ]
> > > 
> > > Second, figure out how to compute it.
> > > 
> > > A quick look at the matlab/octave command
> > > 	plot( exp( linspace(0,1) ) )
> > > shows that the exp function fairly linear over the range [0,1].
> > > It is a good candidate for a lookup table with linear interpolation.
> > > 
> > > With Q.15 ints, then the input range 0:1 is all you need.
> > > 
> > > If your input range can be > 1, you could either
> > > extend the range of your lookup table
> > > or
> > > shift the input down to the range 0:1 and use the identity
> > > e^x = ( e^(x/n) )^n.
> > > 
> > > 
> > > Without more analysis, I'm not sure which would be better.
> > > You must decide for yourself.  Analyze it. Model it. Test it. Report 
> > > your results next week.  (I'm not entirely kidding. I want to know. )
> > > 
> > > Exponentiation becomes increasingly non-linear (and non-representable in 
> > > fixed-point ) at higher input values. That might tend to offset the loss 
> > > of precision due to squaring of the interpolated answer.  If you have 
> > > more than a few integer bits of the Q representation, you may want to 
> > > saturate the output answers at some point to avoid horrible loss of 
> > > precision at lower input values.
> > > 
> > > 
> > > Cheers,
> > > Mark
> > 
> > Thank you very much.
> > 
> > look up table is much faster, because I have a great external memory.
> > But my input is varying in a big range(the results are from 1 to
> > 1E18), so, I think I will analyze the two methods you said. I will
> > paste the result.
> > 
> > peace,
> > 
> > Ken
> 
> A common solution to this problem is to use range reduction followed
> by a polynomial fit. For a discussion, see the book COMPUTER
> APPROXIMATIONS, by J.F. Hart.
> 
> There is free assembler code to do exp() in the TI DSPLIB software
> package found at www.ti.com, document # SPRA480B.
> 
> Regards,
> 
> John

gotcha,a very useful document, thank you very much,

peace,

Ken

kenhusier@yahoo.com (Ken) wrote in message news:<46c498f7.0407060822.6197b9b7@posting.google.com>...
> Mark Borgerding <mark@borgerding.net> wrote in message news:<MmmFc.6272$m91.5379@fe1.columbus.rr.com>...
> > Ken wrote:
> > > 
> > > I have a question, if I need to calculate float point exp()function,
> > I assume you mean fixed point exp(), right?
> > 
> > > what should I do? have to change them back to float, right? if then, I
> > > still can't avoid from a floating point arithmatic.
> > > 
> > > Thanx,
> > > Ken
> > 
> > You left out a VERY important piece of info.  What is your fixed point 
> > representation?  I'm assuming Q15 , since implementing an exp() function 
> > for fixed-point representation with more than one or two integer bits 
> > becomes ridiculous.   exp(127) =~ 1e55
> > 
> > 
> > First, decide how you'll represent the return value.
> > 
> > If you're using Q15 ints, the range is [-1:1).
> > exp(-1) = 1/e =~ .367
> > exp(1) = e = ~ 2.718
> > 
> > Assuming you can overcome a constant multiplier, I suggest returning 
> > exp(x)/e to prevent overflow.
> > 
> > The output range becomes [ exp(-2):1 ] =~ [ .13534... : 1 ]
> > 
> > Second, figure out how to compute it.
> > 
> > A quick look at the matlab/octave command
> > 	plot( exp( linspace(0,1) ) )
> > shows that the exp function fairly linear over the range [0,1].
> > It is a good candidate for a lookup table with linear interpolation.
> > 
> > With Q.15 ints, then the input range 0:1 is all you need.
> > 
> > If your input range can be > 1, you could either
> > extend the range of your lookup table
> > or
> > shift the input down to the range 0:1 and use the identity
> > e^x = ( e^(x/n) )^n.
> > 
> > 
> > Without more analysis, I'm not sure which would be better.
> > You must decide for yourself.  Analyze it. Model it. Test it. Report 
> > your results next week.  (I'm not entirely kidding. I want to know. )
> > 
> > Exponentiation becomes increasingly non-linear (and non-representable in 
> > fixed-point ) at higher input values. That might tend to offset the loss 
> > of precision due to squaring of the interpolated answer.  If you have 
> > more than a few integer bits of the Q representation, you may want to 
> > saturate the output answers at some point to avoid horrible loss of 
> > precision at lower input values.
> > 
> > 
> > Cheers,
> > Mark
> 
> Thank you very much.
> 
> look up table is much faster, because I have a great external memory.
> But my input is varying in a big range(the results are from 1 to
> 1E18), so, I think I will analyze the two methods you said. I will
> paste the result.
> 
> peace,
> 
> Ken

A common solution to this problem is to use range reduction followed
by a polynomial fit. For a discussion, see the book COMPUTER
APPROXIMATIONS, by J.F. Hart.

There is free assembler code to do exp() in the TI DSPLIB software
package found at www.ti.com, document # SPRA480B.

Regards,

John

Mark Borgerding <mark@borgerding.net> wrote in message news:<MmmFc.6272$m91.5379@fe1.columbus.rr.com>...
> Ken wrote:
> > 
> > I have a question, if I need to calculate float point exp()function,
> I assume you mean fixed point exp(), right?
> 
> > what should I do? have to change them back to float, right? if then, I
> > still can't avoid from a floating point arithmatic.
> > 
> > Thanx,
> > Ken
> 
> You left out a VERY important piece of info.  What is your fixed point 
> representation?  I'm assuming Q15 , since implementing an exp() function 
> for fixed-point representation with more than one or two integer bits 
> becomes ridiculous.   exp(127) =~ 1e55
> 
> 
> First, decide how you'll represent the return value.
> 
> If you're using Q15 ints, the range is [-1:1).
> exp(-1) = 1/e =~ .367
> exp(1) = e = ~ 2.718
> 
> Assuming you can overcome a constant multiplier, I suggest returning 
> exp(x)/e to prevent overflow.
> 
> The output range becomes [ exp(-2):1 ] =~ [ .13534... : 1 ]
> 
> Second, figure out how to compute it.
> 
> A quick look at the matlab/octave command
> 	plot( exp( linspace(0,1) ) )
> shows that the exp function fairly linear over the range [0,1].
> It is a good candidate for a lookup table with linear interpolation.
> 
> With Q.15 ints, then the input range 0:1 is all you need.
> 
> If your input range can be > 1, you could either
> extend the range of your lookup table
> or
> shift the input down to the range 0:1 and use the identity
> e^x = ( e^(x/n) )^n.
> 
> 
> Without more analysis, I'm not sure which would be better.
> You must decide for yourself.  Analyze it. Model it. Test it. Report 
> your results next week.  (I'm not entirely kidding. I want to know. )
> 
> Exponentiation becomes increasingly non-linear (and non-representable in 
> fixed-point ) at higher input values. That might tend to offset the loss 
> of precision due to squaring of the interpolated answer.  If you have 
> more than a few integer bits of the Q representation, you may want to 
> saturate the output answers at some point to avoid horrible loss of 
> precision at lower input values.
> 
> 
> Cheers,
> Mark

Thank you very much.

look up table is much faster, because I have a great external memory.
But my input is varying in a big range(the results are from 1 to
1E18), so, I think I will analyze the two methods you said. I will
paste the result.

peace,

Ken

Ken wrote:
> 
> I have a question, if I need to calculate float point exp()function,
I assume you mean fixed point exp(), right?

> what should I do? have to change them back to float, right? if then, I
> still can't avoid from a floating point arithmatic.
> 
> Thanx,
> Ken

You left out a VERY important piece of info.  What is your fixed point 
representation?  I'm assuming Q15 , since implementing an exp() function 
for fixed-point representation with more than one or two integer bits 
becomes ridiculous.   exp(127) =~ 1e55


First, decide how you'll represent the return value.

If you're using Q15 ints, the range is [-1:1).
exp(-1) = 1/e =~ .367
exp(1) = e = ~ 2.718

Assuming you can overcome a constant multiplier, I suggest returning 
exp(x)/e to prevent overflow.

The output range becomes [ exp(-2):1 ] =~ [ .13534... : 1 ]

Second, figure out how to compute it.

A quick look at the matlab/octave command
	plot( exp( linspace(0,1) ) )
shows that the exp function fairly linear over the range [0,1].
It is a good candidate for a lookup table with linear interpolation.

With Q.15 ints, then the input range 0:1 is all you need.

If your input range can be > 1, you could either
extend the range of your lookup table
or
shift the input down to the range 0:1 and use the identity
e^x = ( e^(x/n) )^n.


Without more analysis, I'm not sure which would be better.
You must decide for yourself.  Analyze it. Model it. Test it. Report 
your results next week.  (I'm not entirely kidding. I want to know. )

Exponentiation becomes increasingly non-linear (and non-representable in 
fixed-point ) at higher input values. That might tend to offset the loss 
of precision due to squaring of the interpolated answer.  If you have 
more than a few integer bits of the Q representation, you may want to 
saturate the output answers at some point to avoid horrible loss of 
precision at lower input values.


Cheers,
Mark

axl_fugazi@yahoo.com (Pallav) wrote in message news:<baaddb55.0407010748.209d2360@posting.google.com>...
> > Here are some macros I define in kissfft for fixed point math (Q.15).
> > 
> > 
> > 
> > #define smul(a,b) ( (long)(a)*(b) )
> > #define sround( x )  (short)( ( (x) + (1<<14) ) >>15 )
> > 
> > #define S_MUL(a,b) sround( smul(a,b) )
> > 
> > #define C_MUL(m,a,b) \
> >        do{ (m).r = sround( smul((a).r,(b).r) - smul((a).i,(b).i) ); \
> >            (m).i = sround( smul((a).r,(b).i) + smul((a).i,(b).r) ); while(0)
> > 
> > 
> > 
> > Hope this helps,
> > Mark
> 
> Hi Mark,
> 
> Thank you for your post. Yes, I made a typo in the posted code. Rather
> than a comma it should have been a *. I will try out your macros.
> 
> thanks again
> pallav

hi, mark & pallav,...

I have a question, if I need to calculate float point exp()function,
what should I do? have to change them back to float, right? if then, I
still can't avoid from a floating point arithmatic.

Thanx,
Ken

> Here are some macros I define in kissfft for fixed point math (Q.15).
> 
> 
> 
> #define smul(a,b) ( (long)(a)*(b) )
> #define sround( x )  (short)( ( (x) + (1<<14) ) >>15 )
> 
> #define S_MUL(a,b) sround( smul(a,b) )
> 
> #define C_MUL(m,a,b) \
>        do{ (m).r = sround( smul((a).r,(b).r) - smul((a).i,(b).i) ); \
>            (m).i = sround( smul((a).r,(b).i) + smul((a).i,(b).r) ); while(0)
> 
> 
> 
> Hope this helps,
> Mark

Hi Mark,

Thank you for your post. Yes, I made a typo in the posted code. Rather
than a comma it should have been a *. I will try out your macros.

thanks again
pallav

Pallav wrote:
> hello,
> 
> i'm trying to implement some fixed point arithmetic. the base i've
> chosen is 22.10 format. to multiply two fixed point numbers i say the
> following


Here are some macros I define in kissfft for fixed point math (Q.15).



#define smul(a,b) ( (long)(a)*(b) )
#define sround( x )  (short)( ( (x) + (1<<14) ) >>15 )

#define S_MUL(a,b) sround( smul(a,b) )

#define C_MUL(m,a,b) \
       do{ (m).r = sround( smul((a).r,(b).r) - smul((a).i,(b).i) ); \
           (m).i = sround( smul((a).r,(b).i) + smul((a).i,(b).r) ); while(0)



Hope this helps,
Mark

Pallav wrote:
> hello,
> 
> i'm trying to implement some fixed point arithmetic. the base i've
> chosen is 22.10 format. to multiply two fixed point numbers i say the
> following
> 
> #define FRACTION 10   /* about 3 digits of precision  */
> #define FPMUL(a1, a2)  (((long long int) a1, (long long int) a2) >>
> FRACTION)

This doesn't look right.  You've got a comma where I'd expect to see '*'.

FPMUL(2,5)
will evaluate as
(a,b) >> 10
and so the comma will evaluate a, then b, and result will be
b>>10
... probably not what you want.

> #define FLOAT2FP(a) ((long) a * (1 << FRACTION))
> #define INT2FP(a) (a << FRACTION)
> 
> and i have macros to convert from int to fp, float to fp, vice versa.
> 
> so then I define a fixed point value of 40.6
> #define MAGIC_VALUE FLOAT2FP(40.6)
> 
> then i write some test code 
> 
> long fpn;
> double fpy;
> 
> for (i = 0; ; i++)
> {
>  fpn = INT2FP(i); // conver to fixed point
>  fpn = FMUL(fpn, MAGIC_VALUE); // multiply in fixed point

Wait a minute. This doesn't even match your previous defintion
( FMUL != FPMUL )


If you've typoed this by hand, which one is wrong?


I doubt you posted the code you intended.  I'm afraid I can't help much.

-- Mark

hello,

i'm trying to implement some fixed point arithmetic. the base i've
chosen is 22.10 format. to multiply two fixed point numbers i say the
following

#define FRACTION 10   /* about 3 digits of precision  */
#define FPMUL(a1, a2)  (((long long int) a1, (long long int) a2) >>
FRACTION)
#define FLOAT2FP(a) ((long) a * (1 << FRACTION))
#define INT2FP(a) (a << FRACTION)

and i have macros to convert from int to fp, float to fp, vice versa.

so then I define a fixed point value of 40.6
#define MAGIC_VALUE FLOAT2FP(40.6)

then i write some test code 

long fpn;
double fpy;

for (i = 0; ; i++)
{
 fpn = INT2FP(i); // conver to fixed point
 fpn = FMUL(fpn, MAGIC_VALUE); // multiply in fixed point

 fpy = 40.6 * i; // multiply in double format
 
 // compare values here to find out the error.
}

i'm noticing that as i gets large (say 2500), fpn and fpy differ by as
much as 2 and the difference continues to increase (this is mainly due
to the fractional part). this is becomign a problem. unfortunately, i
can't increase the fractional precision of the fixed point. is there
anything to remedy this? perhaps implementing the FPMUL in a better
way? i don't understand why this is happening given that i have 3
digits of precision?

i also tried doing 40 * i + 0.6 * i, that helps a little but not much.
If i could keep the difference less than 0.5 for all i (within the
limit of the size of an int) i'd be happy.

reason is because i use these values to index a fixed-point
sine/cosine tables and because of this difference (with very large i,
the difference can go upto 20), i'm getting the wrong values during
lookup.

thanks
pallav