[GRASS5] r.mapcalc: x^y and exp(x,y)

Glynn,
(cc grass5)

as you are currently on r.mapcalc:
Do you know the difference between x^y and exp(x,y)?
It appears to be doing the same thing. Do you know what is the difference
between those two?

Thanks,

Markus

Markus Neteler wrote:

as you are currently on r.mapcalc:
Do you know the difference between x^y and exp(x,y)?
It appears to be doing the same thing. Do you know what is the difference
between those two?

Both call pow(), but the "^" operator (implemented by power_x() in
math.c) explicitly checks for the case of a negative number raised to
a non-integral power, and returns NULL (using G_set_d_null_value())

However, in this case, pow() will return NaN (and set errno to EDOM),
so there shouldn't be any difference. However, there is; the two cases
return different NaN values: exp(x,y) returns 0x8000000000000 while
x^y returns 0xfffffffffffff. The former isn't NULL according to
G_is_d_null_value(), while the latter is.

BTW, the current behaviour of r.mapcalc3 matches that of exp() in for
both exp() and the "^" operator when x and y are FCEL/DCELL values.
However, another difference is that the "^" operator in r.mapcalc3
returns CELL if both arguments are CELL, while exp() always returns
DCELL.

I'll fix both r.mapcalc and r.mapcalc3 to ensure that the case of a
negative number raised to a non-integral power is always NULL.
However, I suspect that it might be better to use:

  errno = 0;
  result = pow(x, y);
  if (errno != 0)
    SETNULL_D(&result);

rather than the existing test:

  if (x < 0.0 && y != ceil(y))

The latter can misbehave due to differences in precision.

More generally, this raises the question of whether the existing
implementation of G_is_[cd]_null_value() is correct, i.e. whether it
should treat all NaN values as NULL rather than one specific bit
pattern.

--
Glynn Clements <glynn.clements@virgin.net>