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
fixed-point multiply
Started by ●June 30, 2004
Reply by ●June 30, 20042004-06-30
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 pointWait 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
Reply by ●June 30, 20042004-06-30
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 > followingHere 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
Reply by ●July 1, 20042004-07-01
> 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, > MarkHi 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
Reply by ●July 2, 20042004-07-02
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 > pallavhi, 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
Reply by ●July 2, 20042004-07-02
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, > KenYou 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
Reply by ●July 6, 20042004-07-06
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, > MarkThank 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
Reply by ●July 7, 20042004-07-07
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, > > KenA 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
Reply by ●July 8, 20042004-07-08
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, > > Johngotcha,a very useful document, thank you very much, peace, Ken