Search code examples
cfloating-point

What operations and functions on +0.0 and -0.0 give different arithmetic results?


In C, when ±0.0 is supported, -0.0 or +0.0 assigned to a double typically makes no arithmetic difference. Although they have different bit patterns, they arithmetically compare as equal.

double zp = +0.0;
double zn = -0.0;
printf("0 == memcmp %d\n", 0 == memcmp(&zn, &zp, sizeof zp));// --> 0 == memcmp 0
printf("==          %d\n", zn == zp);                        // --> ==          1

Inspire by a @Pascal Cuoq comment, I am looking for a few more functions in standard C that provide arithmetically different results.

Note: Many functions, like sin(), return +0.0 from f(+0.0) and -0.0 from f(-0.0). But these do not provide different arithmetic results. Also the 2 results should not both be NaN.


Solution

  • There are a few standard operations and functions that form numerically different answers between f(+0.0) and f(-0.0).

    Different rounding modes or other floating point implementations may give different results.

    #include <math.h>
    #include <stdio.h>
    
    double inverse(double x) { return 1/x; }
    
    double atan2m1(double y) { return atan2(y, -1.0); }
    
    double sprintf_d(double x) {
      char buf[20];
      // sprintf(buf, "%+f", x);   Changed to e
      sprintf(buf, "%+e", x);
      return buf[0];  // returns `+` or `-`
    }
    
    double copysign_1(double x) { return copysign(1.0, x); }
    
    double signbit_d(double x) {
      int sign = signbit(x);  // my compile returns 0 or INT_MIN
      return sign;
    }
    
    double pow_m1(double x) { return pow(x, -1.0); }
    
    void zero_test(const char *name, double (*f)(double)) {
      double fzp = (f)(+0.0);
      double fzn = (f)(-0.0);
      int differ = fzp != fzn;
      if (fzp != fzp && fzn != fzn) differ = 0;  // if both NAN
      printf("%-15s  f(+0):%-+15e %s  f(-0):%-+15e\n",
          name, fzp, differ ? "!=" : "==", fzn);
    }
    
    void zero_tests(void) {
      zero_test("1/x",             inverse);
      zero_test("atan2(x,-1)",     atan2m1);
      zero_test("printf(\"%+e\")", sprintf_d);
      zero_test("copysign(x,1)",   copysign_1);
      zero_test("signbit()",       signbit_d);
      zero_test("pow(x,-odd)",     pow_m1);;  // @Pascal Cuoq
      zero_test("tgamma(x)",       tgamma);  // @vinc17 @Pascal Cuoq
    #if __STDC_VERSION__ >= 202310  // C23
      zero_test("rsqrt(x)",        rsqrt);
    #endif
    }
    
    int main(void) {
      zero_tests();
    }
    

    Output:

    1/x              f(+0):+inf            !=  f(-0):-inf           
    atan2(x,-1)      f(+0):+3.141593e+00   !=  f(-0):-3.141593e+00  
    printf("%+e")    f(+0):+4.300000e+01   !=  f(-0):+4.500000e+01  
    copysign(x,1)    f(+0):+1.000000e+00   !=  f(-0):-1.000000e+00  
    signbit()        f(+0):+0.000000e+00   !=  f(-0):+1.000000e+00  
    pow(x,-odd)      f(+0):+inf            !=  f(-0):-inf           
    tgamma(x)        f(+0):+inf            !=  f(-0):-inf      
    

    Notes:
    rsqrt(x), AKA 1/sqrt(x) added to C23. (Untested).

    double zero = +0.0; memcpy(&zero, &x, sizeof x) can show x is a different bit pattern than +0.0 but x could still be a +0.0. I think some FP formats have many bit patterns that are +0.0 and -0.0. TBD.

    This is a self-answer as provided by https://stackoverflow.com/help/self-answer.