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

Question

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.

Katharyn

Note: the memcmp example does not guarantee to find that zn (if we didn't know its value already) is negative zero; it could be an alternative representation of positive zero. – Lyckman

@Matt McNabb Agree about memcmp() as noted in this answer that it lacks definiteness. – Katharyn

Answer 1

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>

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
}

---

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):-2.147484e+09 
pow(x,-odd)      f(+0):+inf             !=  f(-0):-inf           
tgamma(x)        f(+0):+inf             !=  f(-0):+inf  

---

Notes:
tgamma(x) came up == on my gcc 4.8.2 machine, but correctly != on others.

rsqrt(), AKA 1/sqrt() is a maybe future C standard function. May/may not also work.

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.

Bearable

I was pretty impressed when you suggested atan2, and even more when you added sprintf. I haven't been able to think of any further function to distinguish +0 and -0 since then. – Prud

@Pascal Cuoq Your comment inspired this effort. Some hyperbolic functions like csch() may also differentiate, but they do not appear in the C standard. Aside from some systems where 1/0 causes bad things, 1/x seems the simplest test. – Katharyn

pow(0.0, -1) != pow(-0.0, -1) – Prud

tgamma(0.0) != tgamma(-0.0) – Prud

@Pascal Cuoq pow(x,-1), Hmmm: just like 1/x. My gcc 4.8.2 tgamma() returns +INF for both but per this it should return +/-INF. I'll add these 2 later. – Katharyn

@clux cool idea to investigate this further. Thank you for your work. Can you please provide in your answer a little table that maps the different functions to the incarnation of the C standard on which they are introduced (or first provide this difference)? – Kibbutz

@Mark A. Feel open to edit this community wiki answer yourself - or at least start the table. – Katharyn

*(uint64_t *)&fzp != *(uint64_t *)&fzn; is a non completely portable alternative to memcmp. – Becnel

@Becnel That code tests if the bit patterns differ. Floating point could have multiple bit patterns for a +0.0. I have used such double where all otherwise sub-normal codes were treated as 0.0. Certainly your suggest works for binary64, but not double in general. – Katharyn

That's what I meant by non completely portable. It has the same drawbacks as the memcmp alternative plus assumes sizeof(double)==sizeof(uint64_t) and may not work on non IEEE-754 architectures. – Becnel

@Becnel Your point is more clear - an additional reason to not add memcpy() and cast to unsigned compares to the list, which so far are not on this answer's list. – Katharyn

sprintf_d is not protected against buffer overflow. If passed a random value, the 20 character buffer may be too short. You should use snprintf() or use the %+e conversion specifier. Allocating a buffer guaranteed to be large enough for %+f may not be so easy. – Becnel

@Becnel Agree - confident I was originally thinking "%+e" but foolishly coded "%+f". Will fix. – Katharyn

I wonder if 1 / x < 0 is a portable way to test negative float number including negative zero? It works on Windows, not sure if it is portable. – Vallejo

@Vallejo Since C99, signbit() is the best way to " test negative float number including negative zero". 1 / x < 0 is a reasonable alternative, yet not specified by C to get what you want. Note: "Windows" is an OS. The issue is a compiler one, not an OS one. – Katharyn

Answer 2

The IEEE 754-2008 function rsqrt (that will be in the future ISO C standard) returns ±∞ on ±0, which is quite surprising. And tgamma also returns ±∞ on ±0. With MPFR, mpfr_digamma returns the opposite of ±∞ on ±0.

Robin

+1 for tgamma() and the rest. rsqrt(-0.0) returns -INF is surprising. Maybe it is an (over-) simplification of 1.0/sqrt(±0.0) returns 1.0/±0.0? – Katharyn

@Bearable Concerning rSqrt, unfortunately there was no followup to my request during the revision of the old IEEE 754-1985 standard: grouper.ieee.org/groups/754/email/msg03855.html – Robin

Answer 3

I think about this method, but I can't check before weekend, so someone might do some experiments on this, if he/she like, or just tell me that it is nonsense:

• Generate a -0.0f. It should be possible to generate staticly by a assigning a tiny negative constant that underflows float representation.

• Assign this constant to a volatile double and back to float.

  By changing the bit representation 2 times, I assume that the compiler specific standard bit representation for -0.0f is now in the variable. The compiler can't outsmart me there, because a totally other value could be in the volatile variable between those 2 copies.

• compare the input to 0.0f. To detect if we have a 0.0f/-0.0f case

• if it is equal, assign the input volitale double variable, and then back to float.

  I again assume that it has now standard compiler representation for 0.0f

• access the bit patterns by a union and compare them, to decide if it is -0.0f

The code might be something like:

typedef union
{
  float fvalue;
  /* assuming int has at least the same number of bits as float */
  unsigned int bitpat;
} tBitAccess;

float my_signf(float x)
{
  /* assuming double has smaller min and 
     other bit representation than float */

  volatile double refitbits;
  tBitAccess tmp;
  unsigned int pat0, patX;

  if (x < 0.0f) return -1.0f;
  if (x > 0.0f) return 1.0f;

  refitbits = (double) (float) -DBL_MIN;
  tmp.fvalue = (float) refitbits;
  pat0 = tmp.bitpat;

  refitbits = (double) x; 
  tmp.fvalue = (float) refitbits;
  patX = tmp.bitpat;

  return (patX == pat0)? -1.0f : 1.0f;

}

• It is not a standard function, or an operator, but a function that should differentiate between signs of -0.0 and 0.0.
• It based (mainly) on the assumption that the compiler vendor does not use different bit patterns for -0.0f as result of changing of formats, even if the floating point format would allow it, and if this holds, it is independent from the chosen bit pattern.
• For a floating point formats that have exact one pattern for -0.0f this function should safely do the trick without knowledge of the bit ordering in that pattern.
• The other assumptions (about size of the types and so on) can be handled with precompiler switches on the float.h constants.

---

Edit: On a second thought: If we can force the value comparing to (0.0 || -0.0) below the smallest representable denormal (subnormal) floating point number or its negative counterpart, and there is no second pattern for -0.0f (exact) in the FP format, we could drop the casting to volatile double. (But maybe keep the float volatile, to ensure that with deactivated denormals the compiler can't do any fancy trick, to ignore operations, that reduce any further the absolut value of things comparing equal to 0.0.)

The Code then might look like:

typedef union
{
  float fvalue;
  /* assuming int has at least the same number of bits as float */
  unsigned int bitpat;
} tBitAccess;

float my_signf(float x)
{

  volatile tBitAccess tmp;
  unsigned int pat0, patX;

  if (x < 0.0f) return -1.0f;
  if (x > 0.0f) return 1.0f;

  tmp.fvalue = -DBL_MIN;

  /* forcing something compares equal to 0.0f below smallest subnormal 
     - not sure if one abs()-factor is enough */
  tmp.fvalue = tmp.fvalue * fabsf(tmp.fvalue);
  pat0 = tmp.bitpat;

  tmp.fvalue = x; 
  tmp.fvalue = tmp.fvalue * fabsf(tmp.fvalue);
  patX = tmp.bitpat;

  return (patX == pat0)? -1.0f : 1.0f;

}

This might not work with fancy rounding methods, that prevent rounding from negative values towards -0.0.

Kibbutz

I see, code is forcing/coaxing any -0.0f to become the -0.0f. A hole could exist with some embedded processors that use a floating point format like IEEE, but when the exponent was 0, the value was 0.0 regardless of the mantissa (no sub-normal numbers) - thus multiple bit patterns that realize 0.0. This same FP had float == double, but I forget if it had a true -0.0. Still nice idea. +1 – Katharyn

compare the input to 0.0f - positive and negative zero must compare equal (if using a comparison operator) – Lyckman

@Matt McNabb You are correct, the last patX == pat0 needs to be something like memcmp() instead as the logical compare == compares the arithmetic values and not the bit patterns. – Katharyn

well, chux: regarding the compare: patX and pat0 is the content of the storage of the float, re-interpreted as unsigned integer. Given the same length of float and unsigned it, patX and pat0 represent the bitpatterns. But now unsigned integers with different bitpatterns are different arithemetic values. Therefore "==" is correct. @MattMcNabb I "compare" the input to 0.0 to detect, if its (0.0 || -0.0) because any other value has another bit pattern to. In the code this is done by early return. – Kibbutz

@clux To be sure, that the denormals won't be a problem, we could force the value (if it compares to 0.0f) below the smallest subnormal pattern. If I renember right, the difference between any two subnormals is equi-distant, but anyway, it should be easy to generate such a pattern by something like: x = x * fabs(x); or a chain of this. Same for the small Constant that provides the reference pattern. ...= -DBL_MIN * DBL_MIN; - so if there is ONE 0.0f pattern, we can construct it. And then the function is exact. – Kibbutz

Answer 4

Not exactly an answer to the question, but can be useful to know:

Just faced the case a - b = c => b = a - c, which fails to hold if a is 0.0 and b is -0.0. We have 0.0 - (-0.0) = 0.0 => b = 0.0 - 0.0 = 0.0. The sign is lost. The -0.0 is not recovered.

Taken from here.

Garman

Interesting. Agree the sign is lost, yet arithmetic compare b(initial) == b(later) remains true. – Katharyn