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

[Glynn Clements]

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 at virgin.net>