Floating point values with sin() and sinf()

Question

Why do I get the same value from sin(x) and sinf(x) as in (8) and (9) below? Why do I get different values of [ x - sinf(x) ] from two different implementations, (2) and (3) below? Why (6) gives the same result as (5)?

I am using g++ (Ubuntu/Linaro 4.6.3-1ubuntu5) 4.6.3. I used -O0 to disable optimization.

Background: I was tracking down a bug in my program where I am required to use float throughout the program because it will be ported to an embedded system. However, I am currently debugging on Ubuntu on my PC just because it is convenient. I found the operation like (x-s) caused inaccuracy when x was small. That made me think it must be due to loss of significant digits by catastrophic cancellation. However, when I replaced a variable, s, with sinf(x), the problem of inaccuracy did not occur (as seen in (2) vs (3)). I could guess sinf() might be implemented as the same thing as sin(). If so, why does the explicit casting to float have no effect as in (4) and (5). Now I am puzzled.

int main()
{
    unsigned long int xx(0x3d65c2f2);
    float x(*reinterpret_cast<float*>(&xx));
    float s(sinf(x));

    printf("( 1) x                         = %.10e\n", x);
    printf("( 2) x - s                     = %.10e\n", x-s);
    printf("( 3) x - sinf(x)               = %.10e\n", x-sinf(x)); // Why is it different from (2)?
    printf("( 4) x - float(sinf(x))        = %.10e\n", x-float(sinf(x))); // Compare with (3). Why casting has no effect?
    printf("( 5) float(x) - float(sinf(x)) = %.10e\n", float(x)-float(sinf(x))); // Compare with (3). Why casting has no effect?
    printf("( 6) x - sin(x)                = %.10e\n", x - sin(x));

    printf("( 7) s                         = %.10e\n", s);
    printf("( 8) sinf(x)                   = %.10e\n", sinf(x));
    printf("( 9) sin(x)                    = %.10e\n", sin(x)); // Compare with (8). Is sinf() identical to sin()?
    printf("(10) float(sinf(x))            = %.10e\n", float(sinf(x))); // Compare with (8). Why casting has no effect?

    double s_df(sinf(x));
    double s_dd(sin(x));
    float s_fd(sin(x));
    float s_ff(sinf(x));
    printf("(20) s_df               = %.10e\n", s_df);
    printf("(21) s_dd               = %.10e\n", s_dd); // Compare with (20). Is sinf() identical to sin()?
    printf("(22) s_fd               = %.10e\n", s_fd);
    printf("(23) s_ff               = %.10e\n", s_ff);

    return 0;
} 

Here is the output:

$ make && ./main
g++ main.cc -Wall  -c -o main.o -O0
g++ -o main main.o  
( 1) x                         = 5.6094117463e-02
( 2) x - s                     = 2.9411166906e-05
( 3) x - sinf(x)               = 2.9412529899e-05
( 4) x - float(sinf(x))        = 2.9412529899e-05
( 5) float(x) - float(sinf(x)) = 2.9412529899e-05
( 6) x - sin(x)                = 2.9412529899e-05
( 7) s                         = 5.6064706296e-02
( 8) sinf(x)                   = 5.6064704933e-02
( 9) sin(x)                    = 5.6064704933e-02
(10) float(sinf(x))            = 5.6064704933e-02
(20) s_df               = 5.6064704933e-02
(21) s_dd               = 5.6064704933e-02
(22) s_fd               = 5.6064706296e-02
(23) s_ff               = 5.6064706296e-02

Answer 1

In C++, sin has an overload float sin(float f). And overload resolution is done on argument type, rather than return type. To force the use of double sin(double d) you need to cast the argument: sin(static_cast<double>(x)).

(2) vs (3): the FP standard allows implementations to store intermediate results with greater accuracy than the final result. So the value of s need not be exactly the same as the intermediate result for sin(f) in (3).

A lot of this is dependent on your compiler, compiler settings and hardware. For example, if I run your code on my system I get:

( 1) x                         = 5.6094117463e-02
( 2) x - s                     = 2.9411166906e-05
( 3) x - sinf(x)               = 2.9411166906e-05
( 4) x - float(sinf(x))        = 2.9411166906e-05
( 5) float(x) - float(sinf(x)) = 2.9411166906e-05
( 6) x - sin(x)                = 2.9412529899e-05
( 7) s                         = 5.6064706296e-02
( 8) sinf(x)                   = 5.6064706296e-02
( 9) sin(x)                    = 5.6064704933e-02
(10) float(sinf(x))            = 5.6064706296e-02
(20) s_df               = 5.6064706296e-02
(21) s_dd               = 5.6064704933e-02
(22) s_fd               = 5.6064706296e-02
(23) s_ff               = 5.6064706296e-02