4

For logical and integer vectors, there is !anyNA(x). For double vectors, there is length(x) == 0L || all(is.finite(range(x))).

For complex vectors, is there an analogous, O(1)-allocating expression equivalent to all(is.finite(x)), in base R?

I considered is.finite(sum(x)), but that could overflow DBL_MAX and in that case give a false negative result ...

> r <- .Machine$double.xmax
> x <- complex(2L, r, r)
> all(is.finite(x))
[1] TRUE
> is.finite(sum(x))
[1] FALSE

I should add:

> str(.Machine)
List of 19
 $ double.eps           : num 2.22e-16
 $ double.neg.eps       : num 1.11e-16
 $ double.xmin          : num 2.23e-308
 $ double.xmax          : num 1.8e+308
 $ double.base          : int 2
 $ double.digits        : int 53
 $ double.rounding      : int 5
 $ double.guard         : int 0
 $ double.ulp.digits    : int -52
 $ double.neg.ulp.digits: int -53
 $ double.exponent      : int 11
 $ double.min.exp       : int -1022
 $ double.max.exp       : int 1024
 $ integer.max          : int 2147483647
 $ sizeof.long          : int 8
 $ sizeof.longlong      : int 8
 $ sizeof.longdouble    : int 8
 $ sizeof.pointer       : int 8
 $ sizeof.time_t        : int 8
Mikael Jagan
  • 9,012
  • 2
  • 17
  • 48
  • 1
    Can you provide a simple [reproducible example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) that can be used to test and verify possible solutions. – MrFlick Jun 26 '23 at 20:13
  • 1
    This is an interesting question in any case, but can you clarify if there's a practical use case? (I suspect the answer is "no, you'll have to write Rcpp code", but can't rule out someone being clever) – Ben Bolker Jun 26 '23 at 20:53
  • @BenBolker Mainly my curiosity. In practice, if I've already allocated `2 * n * sizeof(double)` bytes, then I probably don't mind allocating `n * sizeof(int)` bytes. On the other hand, I am assuming something about the ratio of `sizeof(double)` and `sizeof(int)` that the standards probably don't guarantee. – Mikael Jagan Jun 26 '23 at 21:09
  • LAPACK doesn't specify how IEEE special values are handled, so sometimes we want to test for them. Testing in R rather than C (where the LAPACK routine in called) is nice for users who want to see where the error is signaled without looking at C code. But to justify testing in R, one would like the R level test to be as efficient as possible. – Mikael Jagan Jun 26 '23 at 21:22
  • 1
    Does IEEE define finite for complex values? Is it just if real and imaginary components are both finite? – qwr Jun 26 '23 at 21:57
  • I'm not sure about O(1), but what about `!is.na(sum(x))` – AndS. Jun 27 '23 at 00:57
  • @qwr AFAIK, IEEE does not specify the binary format of a complex infinity. Hence in standard C we have `a+inf*i` or `inf+b*i` (represented as an array of two floats) rather than some macro constant analogous to `INFINITY`. – Mikael Jagan Jun 27 '23 at 01:33
  • @AndS. Not equivalent to `all(is.finite(x))`, e.g., for `x <- Inf+0i`. – Mikael Jagan Jun 27 '23 at 01:40
  • Anything wrong with `is.finite(range(Re(x))) && is.finite(range(Im(x)))` – Joseph Wood Jun 27 '23 at 21:54
  • 1
    @JosephWood Assuming that you mean `&` and not `&&`: `Re(x)` and `Im(x)` are each O(n) allocations. That expression would be equivalent to but _less_ efficient than `all(is.finite(x))`. – Mikael Jagan Jun 27 '23 at 22:46

2 Answers2

2

We could use mean(x):

is.finite(mean(x))
isTRUE(mean(x)-mean(x)==0)
is.finite(Re(mean(x)))&is.finite(Im(mean(x)))

Testing:

r <- .Machine$double.xmax
x <- complex(2L, r, r)
is.finite(mean(x))
[1] TRUE
isTRUE(mean(x)-mean(x)==0)
[1] TRUE
is.finite(Re(mean(x)))&is.finite(Im(mean(x)))
[1] TRUE

x2 <- c(x,Inf)
is.finite(mean(x2))
[1] FALSE
isTRUE(mean(x2)-mean(x2)==0)
[1] FALSE
is.finite(Re(mean(x2)))&is.finite(Im(mean(x2)))
[1] FALSE
one
  • 3,121
  • 1
  • 4
  • 24
  • D'oh - yes, `mean` does not overflow `double` as easily as `sum`. Still, I think it would be possible to devise an example where `is.finite(mean(x)) != all(is.finite(x))`. I'll have to look at the source code. Meanwhile, I've upvoted. – Mikael Jagan Jun 27 '23 at 15:03
  • Actually, `mean` works by computing the sums of the real and imaginary parts in `long double` then dividing by `n`. But, on my system (macOS, arm64), `long double` and `double` have the same width, so `mean` is not safer than `sum`. – Mikael Jagan Jun 27 '23 at 15:23
  • I had the question on why does summation in ```mean``` not overflow and thanks for answering that! Indeed, the code fails on my machine when ```x <- complex(4094L, r, r)```. – one Jun 27 '23 at 16:04
  • FYI, there is [iterative mean](https://www.heikohoffmann.de/htmlthesis/node134.html) that does not overflow. It is not in base R but could be implemented easily. – one Jun 27 '23 at 16:06
1

Update

For O(1) space complexity, maybe we have to collapse the vector into a single value first, and one possible option is

Reduce(\(p, q) p & is.finite(q), x, init = TRUE)

or as @peter861222 did

is.finite(mean(x))

or with your custom functions for example

f1 <- function(v) {
    if (length(v) == 1) {
        return(is.finite(v))
    }
    is.finite(v[1]) & Recall(v[-1])
}

f2 <- function(v) {
    r <- TRUE
    for (k in v) {
        r <- r & is.finite(k)
    }
    r
}

Previous Solution (NOT O(1) space complexity)

Probably you can try

> all(1 / x != 0, na.rm = TRUE)
[1] TRUE

and with Inf and -Inf in the complex vector, we can see

> (y <- c(x, Inf, -Inf))
[1] 1.797693e+308+1.797693e+308i 1.797693e+308+1.797693e+308i
[3]           Inf+ 0.000000e+00i          -Inf+ 0.000000e+00i

> all(1 / y != 0, na.rm = TRUE)
[1] FALSE
ThomasIsCoding
  • 96,636
  • 9
  • 24
  • 81
  • how does it handle 0+0i – qwr Jun 26 '23 at 21:49
  • @qwr then we can enable `na.rm = TRUE` in `all`. – ThomasIsCoding Jun 26 '23 at 21:52
  • isn't computing `1/x` going to lead to allocating a new complex vector of length N ... ? – Ben Bolker Jun 26 '23 at 23:11
  • This answer fails the "O(1)-allocating" requirement. The division by `x` and the comparison with 0 each require O(n) allocations. – Mikael Jagan Jun 26 '23 at 23:12
  • @BenBolker aha, yes, you are right. I overlooked the condition of `O(1)` from OP's question. – ThomasIsCoding Jun 27 '23 at 08:19
  • @MikaelJagan what about `Reduce` or `mean` in my update? – ThomasIsCoding Jun 27 '23 at 08:30
  • `body(Reduce)` suggests that it is also O(n). `mean.default` is closer to C code and seems to be the closest approximation to an answer. I'm going to play around to see if I can find an example where `is.finite(mean.default(x)) != all(is.finite(x))`. – Mikael Jagan Jun 27 '23 at 14:56
  • @MikaelJagan I am not sure if `O(1)` is really achievable. and perhaps you can borrows the concepts from `Reduce` and define your own function. I updated my solution with two functions for example, but the recursion approach `f1` not recommended for efficiency purpose. – ThomasIsCoding Jun 27 '23 at 18:32
  • Yes, I think that you are right. But the clever people here have surprised me before, so I figured that I should ask anyway. In practice, I will mostly happily continue to use the non-obfuscating `all(is.finite(x))`. – Mikael Jagan Jun 27 '23 at 19:22