3

I noticed some strange behavior in one of my compute shaders where they would return nan unexpectedly. When investigating closer I found pow to be the culprit:

pow(1, inf) == NaN

From a C/C++ standpoint I would expect pow(1, inf) == 1:

pow(+1, exponent) returns 1 for any exponent, even when exponent is NaN

Is GLSL's pow defined for this input?

I only found the usual exceptions for base < 0, and base == 0 && exponent <= 0 in the specification.

Nobody moving away from SE
  • 9,347
  • 3
  • 45
  • 92
  • 2
    `1/0`, huh? division by zero? infinity seems like the proper result and NaN as the power of infinity would be impossible. – vallentin Aug 10 '15 at 09:10
  • @Vallentin: I can only tell you what I am seeing. No compiler errors (or warnings) for 1/0 and the result of 1 from pow. If you have a better way to [create INF inside GLSL](https://stackoverflow.com/questions/10435253/glsl-infinity-constant) tell me about it. – Nobody moving away from SE Aug 10 '15 at 09:13
  • division by zero doesn't need to throw an error as infinity is an accepted result. Thereby x^infinity again equals NaN which would be expected. What I'm saying is, I can't see what the problem is, as the result makes sense. – vallentin Aug 10 '15 at 09:15
  • 1
    @Vallentin From a mathematical standpoint I would expect 1^infinity to be 1, x^infinity for 0 <= x < 1 to be 0, and x^infinity for x > 1 to be infinity. The only place for NaN IMHO is negative numbers to the power of infinity. – Nobody moving away from SE Aug 10 '15 at 09:18
  • True, though it might be because it can't represent the number and thereby gives NaN, why it doesn't give infinity I don't know. But that's actually a good question then. – vallentin Aug 10 '15 at 09:21
  • @Vallentin: I removed the `1/0` part. I made the integer division error. For true INF values (1.0/0.0) the return value is consistently NaN. (It seems `1/0 == -1` in GLSL, so the expression was `pow(1,-1)` which is cleary 1). – Nobody moving away from SE Aug 10 '15 at 09:27
  • 4
    @Nobody From a mathematical perspective, 1^Infinity is not 1. http://mathforum.org/library/drmath/view/56006.html – nos Aug 10 '15 at 11:19
  • @nos: You are right, I totally forgot about the limit from other directions so I should not motivate this with maths. I expected GLSL to behave like [C](http://en.cppreference.com/w/c/numeric/math/pow)/[C++](http://en.cppreference.com/w/cpp/numeric/math/pow) where cppreference.com states: `pow(+1, exponent) returns 1 for any exponent, even when exponent is NaN`. – Nobody moving away from SE Aug 10 '15 at 12:09
  • Technically C++ does not require that at all, it's a [IEE754 pecularity](http://grouper.ieee.org/groups/754/meeting-materials/2001-07-18-c99.pdf). And GLSL [doesn't do IEEE754](http://stackoverflow.com/a/12223539/15416) – MSalters Aug 10 '15 at 12:33
  • @MSalters: Okay, so another region where the C++ standard gives leeway to the compiler writers? As to GLSL: your link does only reiterated what I already said. I took closer look at the specification and it states: `While encodings are logically IEEE 754, operations (addition, multiplication, etc.) are not necessarily performed as required by IEEE 754` and `Infinities and zeros are generated as dictated by IEEE, but subject to the precisions allowed in the following table`. Those two citations are from 4.1.4 and 4.7.1 of GLSLangSpec.4.40 and (IMHO) no further relevant mentions of NaN are made. – Nobody moving away from SE Aug 10 '15 at 13:04
  • Indeed, C++ doesn't require IEEE754 and often it's not (fully) implemented. – MSalters Aug 10 '15 at 13:07
  • @MSalters: Stop destroying my orderly world view ;) Do you have some sources that state where common C++ implementations differ from IEEE? – Nobody moving away from SE Aug 10 '15 at 13:21
  • Truncating denormals to zero is common, because the CPU doesn't support them in **real** hardware. They'll use microcode instead. Compilers optimizing a+b-b to just a. That kind of stuff. – MSalters Aug 10 '15 at 13:35

2 Answers2

4

@Damon's answer was close but copied a mistake from the spec:

pow(x, y) == exp2(x *  log2 (y))

Which has x and y swapped. It should be

pow(x, y) == exp2(y *  log2 (x))

(test e.g pow(2, 3))

Using the correct formula we get

pow(1, inf) == exp2(inf * log2(1)) == exp2(inf * 0) == exp2(nan) == nan
Nobody moving away from SE
  • 9,347
  • 3
  • 45
  • 92
2

Although compute shaders kind of imply that you use a GPU and GLSL version (> 1.30) that does IEEE compliant floating-point math, this is not something you can take for granted, unluckily. So... here be dragons.

Your reasoning from a C++ point of view quoting the C++ reference is reasonable and correct, but it omits a small but important sentence: "If the implementation supports IEEE floating-point arithmetic (IEC 60559),". This small sentence highlights that it's neither a C++ requirement (just a coincidence since IEEE 754 asks for it, and usually C++ implementations comply with IEEE), and it's not a strict guarantee either (it is perfectly allowable for a C++ implementation to not implement IEEE 754).

Further, GLSL is not C++. This may come as a surprise with arithmetic because although it obviously isn't the same language, one can expect that math still works the same, can't one?

Well, no. pow(x, y) is defined (although 8.2 on which the online man page is based does not explicitly say so) as exp2(x*log2(x)) in GLSL. It inherits its precision from the forementioned expression as stated in 4.7.1. It is not stated that pow(1, x) equals 1.

Now, expanding the above expression with the presumably correct intermediate results, the correct result should be +INF, since log2(+INF) = +INF, and 1*+INF is still +INF, and exp2(+INF) = +INF.
But, the standard doesn't really define that. It only explicitly (un)defines zero exponents.

Damon
  • 67,688
  • 20
  • 135
  • 185