24

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.

chux - Reinstate Monica
  • 143,097
  • 13
  • 135
  • 256
  • 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. – M.M Dec 01 '14 at 03:36
  • @Matt McNabb Agree about `memcmp()` as noted in this [answer](http://stackoverflow.com/a/25332135/2410359) that it lacks definiteness. – chux - Reinstate Monica Dec 01 '14 at 04:49

4 Answers4

22

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.

Community
  • 1
  • 1
chux - Reinstate Monica
  • 143,097
  • 13
  • 135
  • 256
  • 4
    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. – Pascal Cuoq Aug 15 '14 at 20:37
  • @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](http://stackoverflow.com/questions/25310051/safe-floating-point-division/25312364) where `1/0` causes _bad_ things, `1/x` seems the simplest test. – chux - Reinstate Monica Aug 15 '14 at 21:27
  • `pow(0.0, -1) != pow(-0.0, -1)` – Pascal Cuoq Aug 15 '14 at 21:51
  • `tgamma(0.0) != tgamma(-0.0)` – Pascal Cuoq Aug 15 '14 at 21:54
  • @Pascal Cuoq `pow(x,-1)`, Hmmm: just like `1/x`. My gcc 4.8.2 `tgamma()` returns +INF for both but per [this](http://en.cppreference.com/w/c/numeric/math/tgamma) it should return +/-INF. I'll add these 2 later. – chux - Reinstate Monica Aug 15 '14 at 22:06
  • @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)? – Mark A. Aug 19 '14 at 18:27
  • @Mark A. Feel open to edit this community wiki answer yourself - or at least start the table. – chux - Reinstate Monica Aug 19 '14 at 20:43
  • `*(uint64_t *)&fzp != *(uint64_t *)&fzn;` is a non completely portable alternative to `memcmp`. – chqrlie May 14 '15 at 21:29
  • @chqrlie 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](http://en.wikipedia.org/wiki/Double-precision_floating-point_format), but not `double` in general. – chux - Reinstate Monica May 14 '15 at 21:38
  • 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. – chqrlie May 14 '15 at 21:40
  • @chqrlie 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. – chux - Reinstate Monica May 14 '15 at 21:53
  • `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. – chqrlie May 14 '15 at 21:59
  • @chqrlie Agree - confident I was originally thinking `"%+e"` but foolishly coded `"%+f"`. Will fix. – chux - Reinstate Monica May 14 '15 at 22:12
  • 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. – raymai97 Dec 08 '18 at 11:20
  • @raymai97 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. – chux - Reinstate Monica Dec 08 '18 at 15:39
7

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.

vinc17
  • 2,829
  • 17
  • 23
  • +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`? – chux - Reinstate Monica Aug 15 '14 at 22:14
  • 2
    @chux Concerning rSqrt, unfortunately there was no followup to my request during the revision of the old IEEE 754-1985 standard: http://grouper.ieee.org/groups/754/email/msg03855.html – vinc17 Aug 15 '14 at 22:24
1

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.

Mark A.
  • 579
  • 4
  • 13
  • 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 – chux - Reinstate Monica Aug 19 '14 at 19:58
  • `compare the input to 0.0f` - positive and negative zero must compare equal (if using a comparison operator) – M.M Aug 19 '14 at 20:07
  • @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. – chux - Reinstate Monica Aug 19 '14 at 20:45
  • 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. – Mark A. Aug 20 '14 at 03:33
  • @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. – Mark A. Aug 20 '14 at 04:00
1

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.

pmor
  • 5,392
  • 4
  • 17
  • 36