3

I am trying to calculate arkus sin() without math library. It works, however produces different results on different platforms.

Windows and Mac (correct): 5.2359877560e-001 Linux: 5.1532736669e-01

Where is problem?

double my_asin(double k)
{
    double x = k;
    double a = 1;
    double s = x;
    const  double hodnota = 0.00000000001f;
    double s1 = 0;
    double d = x;
    double e = 1.0;
    double y;
    y = x * x;


    while(my_abs(s1 - s) >= hodnota) {
        s1 = s;

        d = d * y;

        e = e * ((a * a) / ((++a) * (++a)) );
        s = s + (d * e);
    }

    return s;

}
tshepang
  • 12,111
  • 21
  • 91
  • 136
user3071818
  • 39
  • 3
  • 2
    You should probably tag the right language, also could be useful to know if each OS are 32 or 64bits. – Alexandre Lavoie Dec 06 '13 at 08:21
  • Also, please post a [SSCCE](http://www.sscce.org/). We need the whole program to understand what is happening. Also, please note the exact platforms you ran this on (OS version, HW platform, processor, compiler, compiler version etc...). – sleske Dec 06 '13 at 08:23
  • And BTW, this is most probably a problem with floating point precision. Read up on [floating point](https://en.wikipedia.org/wiki/Floating_point). Note that there is a whole branch of science dedicated to these problems, [Numerical Analysis](https://en.wikipedia.org/wiki/Numerical_analysis). – sleske Dec 06 '13 at 08:26
  • This is yields undefined behavior the following thread explains why. http://stackoverflow.com/questions/4176328/undefined-behavior-and-sequence-points – hetepeperfan Dec 06 '13 at 08:57

2 Answers2

9

Instruction e = e * ((a * a) / ((++a) * (++a)) ); can produce different results because of undefined behavior

You need to change your code this way:

e *= a * a / ((a + 1) * (a + 2));
a += 2;
Community
  • 1
  • 1
vzhilin
  • 313
  • 2
  • 12
2

This line has undefined effect. 

e = e * ((a * a) / ((++a) * (++a)) );

You need to move one, and ideally both the incrmements to a different line. eg: 

e = e * ((a * a) / ((a+2) * (a+1)) );
a+=2

Though you will need to play with the replacement line until it actually does what you hope.

Btw, this code could easily vary across version of compilers, not just compiler brands and OS's.

RichardPlunkett
  • 2,998
  • 14
  • 14