0

EDIT

SOLVED

Solution was to use the long double versions of sin & cos: sinl & cosl.

It is my first post here, so bear with me :).

I come today here to ask for your input on a small problem that I am having with a C application at work. Basically, I am computing an Extended Kalman Filter and one of my formulas (that I store in a variable) has multiple computations of sin and cos, at least 16 in total in the same line. I want to decrease the time it takes for the computation to be done, so the idea is to compute each cos and sin separately, store them in a variable, and then replace the variables back in the formula. So I did this:

const ComputationType sin_Roll = compute_sin((ComputationType)(Roll));
const ComputationType sin_Pitch = compute_sin((ComputationType)(Pitch));

const ComputationType cos_Pitch = compute_cos((ComputationType)(Pitch));
const ComputationType cos_Roll = compute_cos((ComputationType)(Roll));

Where ComputationType is a macro definition (renaming) of the type Double. I know it looks ugly, a lot of maybe unnecessary castings, but this code is generated in Python and it was specifically designed so....

Also, compute_cos and compute_sin are defined as such:

#define compute_sin(a) sinf(a)
#define compute_cos(a) cosf(a)

My problem is the value I get from the "optimized" formula is different from the value of the original one.

I will post the code of both and I apologise in advance because it is very ugly and hard to follow but the main points where cos and sin have been replaced can be seen. This is my task, to clean it up and optimize it, but I am doing it step by step to make sure I don't introduce a bug.

So, the new value is:

ComputationType newValue = (ComputationType)(((((((ComputationType)-1.0))*(sin_Pitch))+((DT)*((((Dg_y)+((((ComputationType)-1.0))*(Gy)))*(cos_Pitch)*(cos_Roll))+(((Gz)+((((ComputationType)-1.0))*(Dg_z)))*(cos_Pitch)*(sin_Roll)))))*(cos_Pitch)*(cos_Roll))+((((DT)*((((Dg_y)+((((ComputationType)-1.0))*(Gy)))*(cos_Roll)*(sin_Pitch))+(((Gz)+((((ComputationType)-1.0))*(Dg_z)))*(sin_Pitch)*(sin_Roll))))+(cos_Pitch))*(cos_Roll)*(sin_Pitch))+((((ComputationType)-1.0))*(DT)*((((Gz)+((((ComputationType)-1.0))*(Dg_z)))*(cos_Roll))+((((ComputationType)-1.0))*((Dg_y)+((((ComputationType)-1.0))*(Gy)))*(sin_Roll)))*(sin_Roll)));

And the original is:

ComputationType originalValue = (ComputationType)(((((((ComputationType)-1.0))*(compute_sin((ComputationType)(Pitch))))+((DT)*((((Dg_y)+((((ComputationType)-1.0))*(Gy)))*(compute_cos((ComputationType)(Pitch)))*(compute_cos((ComputationType)(Roll))))+(((Gz)+((((ComputationType)-1.0))*(Dg_z)))*(compute_cos((ComputationType)(Pitch)))*(compute_sin((ComputationType)(Roll)))))))*(compute_cos((ComputationType)(Pitch)))*(compute_cos((ComputationType)(Roll))))+((((DT)*((((Dg_y)+((((ComputationType)-1.0))*(Gy)))*(compute_cos((ComputationType)(Roll)))*(compute_sin((ComputationType)(Pitch))))+(((Gz)+((((ComputationType)-1.0))*(Dg_z)))*(compute_sin((ComputationType)(Pitch)))*(compute_sin((ComputationType)(Roll))))))+(compute_cos((ComputationType)(Pitch))))*(compute_cos((ComputationType)(Roll)))*(compute_sin((ComputationType)(Pitch))))+((((ComputationType)-1.0))*(DT)*((((Gz)+((((ComputationType)-1.0))*(Dg_z)))*(compute_cos((ComputationType)(Roll))))+((((ComputationType)-1.0))*((Dg_y)+((((ComputationType)-1.0))*(Gy)))*(compute_sin((ComputationType)(Roll)))))*(compute_sin((ComputationType)(Roll)))));

What I want is to get the same value as in the original formula. To compare them I use memcmp.

Any help is welcome. I kindly thank you in advance :).

EDIT

I will post also the values that I get.

New value : -1.2214615708217025e-005

Original value : -1.2214615708215651e-005

They are similar up to a point, but the application is safety critical and it is necessary to validate the results.

Eduard Palko Mate
  • 133
  • 7
  • 1
    How different are the results? Are they completely different or is there just a rounding error? – Blaze Nov 30 '18 at 07:41
  • `precision up to 10^-38` a `double` variable has a precision of approx 15 decimal places... https://stackoverflow.com/questions/9999221/double-precision-decimal-places – Rishikesh Raje Nov 30 '18 at 07:49
  • So you are explicitly using the float version of sin (sinf). Floats have about 5-6 digits of precision, I believe, so what does the expected precision of 10^-38 even mean? – Johnny Johansson Nov 30 '18 at 07:51
  • Also, the precision is for one operation. You are getting 13 digits correct here for many operations. – Rishikesh Raje Nov 30 '18 at 07:52
  • I think you need to use `sinl & cosl` as per https://en.cppreference.com/w/c/numeric/math/sin – ntshetty Nov 30 '18 at 07:58
  • Yes, it is quite ambiguous for me too with the 10^-38... I was told that at some point there would be such values. Anyway, the sinl & cosl did the trick. The values are the same! Thank you for your help. – Eduard Palko Mate Nov 30 '18 at 08:01
  • @EduardPalkoMate Mark it answered, it may help someone – ntshetty Nov 30 '18 at 08:07
  • It is likely that sinl and cosl take a lot more time to compute than sinf and cosf. If the main purpose is to reduce computation time, you might end up doing the opposite. – Johnny Johansson Nov 30 '18 at 08:14
  • They would, but I have 30 lines of formulas with sin and cos. By computing them only once, I will still gain time. – Eduard Palko Mate Nov 30 '18 at 08:19
  • `sinl(x)` is for `long double`. Since you say you're using `double`, just use `sin(x)`. – Giovanni Cerretani Nov 30 '18 at 08:28
  • Follow up on the expected precision. Pitch and Roll data over ARINC bus from the IRS will not be that accurate. It is only 6.7E-4 degrees precision. The rest of the system probably won't even use 64bit floating point (graphics system will use 32bit floats or maybe fixed point integers when Pitch and Roll is displayed to pilot or it'll be rounded again when transmitted over data bus). So getting 13 digits correct seems really good. https://aviation.stackexchange.com/questions/31160/is-there-a-such-thing-as-an-arinc-label-for-the-yaw-angle – fdk1342 Nov 30 '18 at 16:41

1 Answers1

0

You can not meet your expectation for a couple of reasons. By altering the code you adjust the machine instructions being used in subtle ways that will impact the final value. For instance if originally it was using fused multiplies and adds and this is no longer happening it will change the result. You don't mention the target architecture. Some architectures retain more than 64bits in the floating point pipeline. These extra bits get rounded when forced into 64bit memory. Again altering how this works will have minor effects on the final output.

fdk1342
  • 2,681
  • 1
  • 12
  • 15