What to do about infinity?

[nshaffer]

Sometimes an algorithm needs a way to refer to positive or negative infinity (e.g., integration bounds) or test if a number is infinite. How are we going to do that? I see three directions we could take

1. Only implement such algorithms for IEEE-compliant real kinds. Use ieee_arithmetic to set and test for infinities.
2. Define special constants inf_{kind} which is IEEE positive infinity for supported kinds and is a special value (maybe huge(x) + 1.0) for non-IEEE kinds. Likewise, testing for infinities either uses ieee_arithmetic or tests for equality with the defined magic number.
3. Create a derived type that represents infinity and behaves as much like an IEEE infinity as possible when used alongside reals.

Option 1 is simple but means that some stdlib functionality just can't be implemented for all a compiler's real kinds.

Option 2 is simple for IEEE-compliant kinds but introduces a lot of undefined behavior for non-IEEE kinds. (e.g.: If k is a non-IEEE real kind, we can't guarantee that inf_k - 1.0_k does the "right" thing semantically.)

Option 3 is complicated and might be a pain to use in practice. However, it's the only way I know of to get close to IEEE infinity semantics with non-IEEE reals.

I'd like to hear opinions and alternatives, if people have them.

[certik]

Can you give an example when this is needed?

[fiolj]

Looking at, for instance scipy source code, I found about a thousand references to np.inf, but I think most (if not all) are to signal the infinite value as stated by @nshaffer. I agree that is very convenient to signal quadrature limits, or step sizes, or similar. But in this case we only need start implementing the constant (for instance as in option 2) and a routine signaling if it is +inf, -inf, or a finite number.
Something like:

  function is_inf(x) result(y)
    implicit none
    integer :: y     ! Returns 0 if x is finite, -1 if is minus-infinite and +1 if is plus-infinite
    ....
  end function is_inf

Of course, we may just want a True/False result, or something different but with a minimal implementation we could keep advancing and later make additions/modifications as more use cases make it necessary.

[certik]

Aren't the infinities disabled in -ffast-math mode? If so, then this might not work.

[fiolj]

@certik you are right, at least gfortran (silently) disables the ieee infinities, in the sense that it compiles and run but, for instance, ieee_is_finite(x) will give the wrong result.

Additionally, the question on when we need an infinite value could gives a good frame to decide if/what to implement.

As said before, scipy has about 1000 references to np.inf. From these only about 400 are in the code (not in tests, or documentation) and, as far as I can tell it is mainly used to initialize a value that will be tested or signal that a result is infinite (if (s == 0) logs = -inf)
In Scipy it is used in:

• Integration routines
• Optimization routines
• stats (mainly continuum distributions)
  Some uses include definition of domains or regions, initialization of steps, or other values that will be updated

In GSL-2.5 there are about 50 references to GSL_POSINF and they are mostly used in continuum distributions. They use when they can C99X INFINITY, when not the equivalent of huge(x), and fallback to ieee by computing the value of 1.0/0.0.

With all that, I am not sure that we need to define an infinite. I've found it handy for integration domains, but my implementation was rather simple and based in huge(x).

[nshaffer]

@certik Defining quadrature domains was my motivating use-case for this issue. Another use would be constraining curve fitting parameters, e.g., SciPy's scipy.optimize.curve_fit. For these cases, there are simple workarounds to not having a proper, reliable infinity. Sometimes huge(x) is a good enough replacement. Sometimes using character variables or logicals as flags is good enough. Sometimes you write separate routines with different parameter lists to handle cases when it makes sense for an argument to be infinite. From both the library author's point of view and the user's point of view, all these solutions look different and are maintained differently.

For now, maybe it makes sense to let individuals do whatever seems natural for the features they're implementing. Eventually, we will probably want to adopt a uniform approach, though.