CodeWarrior for MCU中的主题 Re: cw v6.2 math precision

cw v6.2 math precision

Sanahuja — Thu, 07 May 2009 04:05:48 GMT

Hello,

I'm working with cw 6.2, and trying do make some trigonometric calculations. The target is a MCF51qe.

For exemple, when I do this:

sin((double)0.86217507327430);

I obtain 0.75926020631887015.

In matlab, the same sinus gives 0.75925986842193.

The difference is not so big, but I have difference for each sin/cos calculation, and the result is used to found the distance between 2 GPS points. The distance is given using an acos function but the argument of the acos is >1, due to precision problems; so it can't be calculated. I'm sure of the formulas I'm using, as I checked with a computer and the problem doesn't appear (argument of acos is always<1).

In cw, I have included this libs:

fp_coldfire_nodiv.a

C_4i_CF_RegABI_MATH_MSL.a

C_4i_CF_RegABI_Runtime.a

I've tried other libs but without success.

Someone has an idea?

Thanks for your help!

Re: cw v6.2 math precision

kef — Thu, 07 May 2009 20:02:41 GMT

0.5e-6 relative error is huge, not "not so big". That's normal performance for single precision FP, not for double!

Re: cw v6.2 math precision

Sanahuja — Thu, 04 Jun 2009 16:39:13 GMT

Nobody has an idea?

If I compute the sinus in single precision, there is always a different between mcu and matlab result.

Thanks

Re: cw v6.2 math precision

kef — Thu, 04 Jun 2009 19:45:06 GMT

Something is wrong with CW library, you should do a service request. If you can't wait for bugfixes, then you may try using third party math functions. Something like cephes lib at www.netlib.org should be more precise.

Re: cw v6.2 math precision

Sanahuja — Fri, 05 Jun 2009 13:46:15 GMT

Thanks for your answer Kef. I will take a look at netlib.org and do a service request.

Re: cw v6.2 math precision

Lundin — Fri, 05 Jun 2009 19:30:49 GMT

I don't think the C standard enforces double to be 64 bytes. When creating a CW project you can decide whether you want double to be 32 bit (same as float) or 64 bit.

Maybe this is the case here?

Re: cw v6.2 math precision

Sanahuja — Fri, 05 Jun 2009 20:36:06 GMT

When I have created the project, to the question:

"select if you want floating point support" I choose "float is ieee64, double is ieee64" .

this includes the foloowing libraries:

C_4i_CF_RegABI_MATH_MSL.a, C++_4i_CF_RegABI_MATH_MSL.a, Cpp_4i_CF_RegABI_Runtime.a

If I choose "no" to the question, the following libraries are included:

C_4i_CF_RegABI_MSL.a, C++_4i_CF_RegABI_MSL.a, Cpp_4i_CF_RegABI_Runtime.a

Re: cw v6.2 math precision

Emac — Thu, 29 Oct 2020 09:28:48 GMT

I can see I neglected to handle to negative (error case)

  if(x < 0.0)  {    return(MY_DOUBLE_ERROR_IS_INFINITY);  }

Or you can use identities to make x = x + 2pi if x < 0

Re: cw v6.2 math precision

Sanahuja — Thu, 18 Jun 2009 16:20:41 GMT

Thanks for your answer Emac.

I will try your sinus function. I also have to check if it is fast enough, as I have to do real time caclulation of GPS position.

Re: cw v6.2 math precision

Emac — Fri, 19 Jun 2009 02:10:09 GMT

If need be, the number of terms can easily be scaled back if you need faster operation at the cost of less accuracy. The current number of terms matched your matlab calculated value but perhaps fewer terms would work as well. Also the constants could be better defined. The source of my numbers and constants was Excell - not known as being the ultimate mathmatical standard.

Good luck!

Re: cw v6.2 math precision

Emac — Thu, 29 Oct 2020 09:28:50 GMT

Final code

/*This formula uses taylor integration expansion and the identity function to evaluatesine over the range of 0 to pi (and by trig identities, -2pi to 2 pi)error should be less than ~ 5e-10*/double mysin(double x){  double identityMatrix[]={-1.00000000000000000E+00,+1.00000000000000000E+00,-1.00000000000000000E+00,+1.00000000000000000E+00,-1.00000000000000000E+00,+1.00000000000000000E+00,-1.00000000000000000E+00};  double multiplierMatrix[] ={5.49450549450549450E-03, //  1/(14*13)7.57575757575757575E-03, //  1/(12*11)1.11111111111111111E-02, //  1/(10*9)1.78571428571428571E-02, //  1/(8*7)3.33333333333333333E-02, //  1/(6*5)8.33333333333333333E-02, //  1/(4*3)5.00000000000000000E-01  //  1/(2*1)};  double result;  double finalSign;  double innerTerm;  const double pi_over_2 = 1.57079632679490000E+00;  const double pi        = 3.14159265358979000E+00;  const double twoPi     = 6.28318530717959000E+00;#define MY_DOUBLE_ERROR_IS_INFINITY 0x7ff0000000000000  //handle negative numbers  if(x < 0.0)  {    x= x + twoPi;  }  //how do you want to handle an error (throw exception?)  if(x >= twoPi)  {    return(MY_DOUBLE_ERROR_IS_INFINITY);  }  if(x < 0.0)  {    return(MY_DOUBLE_ERROR_IS_INFINITY);  }  if( x < pi)  {    innerTerm = x -  pi_over_2;    innerTerm = innerTerm*innerTerm;    finalSign = (double) -1.0;  }else  {    innerTerm = x - pi -  pi_over_2;    innerTerm = innerTerm*innerTerm;    finalSign = (double) +1.0;  }  result = multiplierMatrix[0]*innerTerm + identityMatrix[0];  result = multiplierMatrix[1]*innerTerm*result + identityMatrix[1];  result = multiplierMatrix[2]*innerTerm*result + identityMatrix[2];  result = multiplierMatrix[3]*innerTerm*result + identityMatrix[3];  result = multiplierMatrix[4]*innerTerm*result + identityMatrix[4];  result = multiplierMatrix[5]*innerTerm*result + identityMatrix[5];  result = multiplierMatrix[6]*innerTerm*result + identityMatrix[6];  result = result * finalSign;  return(result);}

I don't know if moderator could link to this code instead of broken one from before

Re: cw v6.2 math precision

Sanahuja — Mon, 22 Jun 2009 17:14:50 GMT

Thanks again Emac.

Your function is working well; I will search for the other trigonometric functions(cos, tan,acos,asin,atan), as I have the same problem with these functions.

I don't had the time to check the speed, but I will do it and let you know.

Re: cw v6.2 math precision

Sanahuja — Thu, 29 Oct 2020 09:28:52 GMT

Here is the response of the service request:

Dear Guillaume,

You are right, the MCF51QE128 isn't as accurate like the matlab
software, the reason is very simple, our mcu doesn't have a
multiplication or division unit, therefore there's the need to create
software routines (libraries) that implements such instructions.

That's the reason why matlab is more accurate, because at matlab all the
processing of information occurs at your computer, the MCF51QE isn't as
powerful like your computer core, but is cheaper and easier to
understand.

For them, the problem is not the library but the MCF51...

Re: cw v6.2 math precision

kef — Mon, 22 Jun 2009 17:49:41 GMT

Then they are not telling the thruth. MCF51 is as accurate as any other computer in the world. Inaccurate is their sin() routine. Their double precision software mul, add, div etc are fine (I checked that if it matters).

It is not true also that MCF51 doesn't have hardware multiplication. It has. It doesn't have division, but division emulation on Coldfire V1 is pretty fast. Of course slower than with hardware division, but not zilions of times.

I think you need accurate results first, then the execution speed. If you would need speed, you could use single precision of maybe fixed point routines.

Re: cw v6.2 math precision

Sanahuja — Mon, 22 Jun 2009 19:17:57 GMT

I agree with you Kef.

I have re-opened the service request.

Re: cw v6.2 math precision

Emac — Thu, 29 Oct 2020 09:28:54 GMT

In case you need a cos algorithem in the mean time.

/*This formula uses taylor integration expansion and the identity function to evaluatecosine over the range of pi/2 to 3pi/2 (and by trig identities, -2pi to 2 pi)error should be less than ~ 1e-10Inspiration from pg 81 tailor series expansion from Numerical Methods fo Engineersby Chapra & Canale Mcgraw Hill publishersa = pi (in order to use trig identities)*/double mycos(double x){  const double multiplierMatrix[] ={-5.49450549450549450E-03, //  -1/(14*13)-7.57575757575757576E-03, //  -1/(12*11)-1.11111111111111111E-02, //  -1/(10*9)-1.78571428571428571E-02, //  -1/(8*7)-3.33333333333333333E-02, //  -1/(6*5)-8.33333333333333333E-02, //  -1/(4*3)-5.00000000000000000E-01 //  -1/(2*1)};  double result;  double finalSign;  double innerTerm;  const double pi_over_2       = 1.57079632679490000E+00;  const double pi              = 3.14159265358979000E+00;  const double three_pi_over_2 = 4.71238898038469000E+00;  const double twoPi           = 6.28318530717959000E+00;#define MY_DOUBLE_ERROR_IS_INFINITY 0x7ff0000000000000  //handle negative numbers  if(x < 0.0)  {    x= x + twoPi;  }  //how do you want to handle an error (throw exception?)  if(x >= twoPi)  {    return(MY_DOUBLE_ERROR_IS_INFINITY);  }  if(x < 0.0)  {    return(MY_DOUBLE_ERROR_IS_INFINITY);  }  if( x < pi_over_2  )  {    innerTerm = x*x;    finalSign = (double) +1.0;  }else if( x < three_pi_over_2 )  {    innerTerm = x - pi;    innerTerm = innerTerm*innerTerm;    finalSign = (double) -1.0;  }else  {    innerTerm = x - twoPi;    innerTerm = innerTerm*innerTerm;    finalSign = (double) +1.0;  }  result = multiplierMatrix[0]*innerTerm        + 1.0;  result = multiplierMatrix[1]*innerTerm*result + 1.0;  result = multiplierMatrix[2]*innerTerm*result + 1.0;  result = multiplierMatrix[3]*innerTerm*result + 1.0;  result = multiplierMatrix[4]*innerTerm*result + 1.0;  result = multiplierMatrix[5]*innerTerm*result + 1.0;  result = multiplierMatrix[6]*innerTerm*result + 1.0;  result = result * finalSign;  return(result);}

Re: cw v6.2 math precision

Emac — Thu, 25 Jun 2009 00:01:14 GMT

I don't know about your compiler library (I use 9s12 ) but the same people most likely did yours as did mine,

I found a conditional compilation statement the adds in additional precision

#if !_SIN_HUGE

...

#else

#ifdef EXTENDED_PRECISION_COSSIN

... //some additional precision here

#endif

... //more precision than original here

#endif

add this to your compiler options or define it in you first file and maybe the precision will come in.

#define _SIN_HUGE

#define EXTENDED_PRECISION_COSSIN

note you may get enough accuracy from _SIN_HUGE alone and you may not need the EXTENDED_PRECISION_COSSIN.

To me the "_" in "_SIN_HUGE" means that it is a system generated definition so you may have to check you options instead of forcing it like I did above.

Check your installation and look at the source file for the math.c file and dig in it to see if you have a conditional compile as well. Then be sure to do a clean compile by removing all of the object files first.

... (my path for example) ...

C:\PROGRA~1\FREESC~1\CWFORH~1.5\lib\hc12c\src

note math.c is for double

mathf.c is for float

Re: cw v6.2 math precision

Sanahuja — Mon, 29 Jun 2009 14:12:28 GMT

Thanks for your cosinus finction Emac, it will be usefull.

I'm quite lost with the libraries, as I'm not used to work with this; but I'm still looking if there is something about precision.

Thanks for your help!

Re: cw v6.2 math precision

ShivaSparr — Sat, 14 Aug 2010 22:20:00 GMT

Hi Emac,

Thanks for posting the codes of sin, cos functions(double precision). Could u please post the code of acos function(double precision).

Waiting for ur reply.

Thanks in advance

Re: cw v6.2 math precision

ShivaSparr — Sat, 14 Aug 2010 22:25:57 GMT

Hi Sanahuja,

I am also facing the similar problem which u were facing, i want to find the distance between 2 GPS points.

I have used the code of sin, cos functions(double precision),which Emac had posted & its working. But i also need the acos function(double precision). So please post the code acos function.

Waiting for ur reply.

Thanks in advance