frollvar and frollsd (#7361)

* frollvar and frollsd

* export

* news entry for var and sd

* codecov fix

* handle var < 0 spotted by Ben

* elaborate more on fall back to slow when NAs are in the input

* align to changes proposed in #7376

* reuse existing macro as per Ben suggestion

Author: jangorecki

NAMESPACE

@@ -58,6 +58,8 @@ export(frollmax)
 export(frollmin)
 export(frollprod)
 export(frollmedian)
+export(frollvar)
+export(frollsd)
 export(frollapply)
 export(frolladapt)
 export(nafill)

NEWS.md

@@ -246,7 +246,7 @@
     #9: 2025-09-22     9         8           9.0
     ```
 
-19. New rolling functions: `frollmin`, `frollprod` and `frollmedian`, have been implemented, towards [#2778](https://github.com/Rdatatable/data.table/issues/2778). Thanks to @jangorecki for implementation. Implementation of rolling median is based on a novel algorithm "sort-median" described by [@suomela](https://github.com/suomela) in his 2014 paper [Median Filtering is Equivalent to Sorting](https://arxiv.org/abs/1406.1717). "sort-median" scales very well, not only for size of input vector but also for size of rolling window.
+19. Other new rolling functions: `frollmin`, `frollprod`, `frollmedian`, `frollvar` and `frollsd`, have been implemented, resolving long standing issue [#2778](https://github.com/Rdatatable/data.table/issues/2778). Thanks to @jangorecki for implementation. Implementation of rolling median is based on a novel algorithm "sort-median" described by [@suomela](https://github.com/suomela) in his 2014 paper [Median Filtering is Equivalent to Sorting](https://arxiv.org/abs/1406.1717). "sort-median" scales very well, not only for size of input vector but also for size of rolling window.
     ```r
     rollmedian = function(x, n) {
       ans = rep(NA_real_, nx<-length(x))

R/froll.R

@@ -216,3 +216,9 @@ frollprod = function(x, n, fill=NA, algo=c("fast","exact"), align=c("right","lef
 frollmedian = function(x, n, fill=NA, algo=c("fast","exact"), align=c("right","left","center"), na.rm=FALSE, has.nf=NA, adaptive=FALSE, partial=FALSE, give.names=FALSE, hasNA) {
   froll(fun="median", x=x, n=n, fill=fill, algo=algo, align=align, na.rm=na.rm, has.nf=has.nf, adaptive=adaptive, partial=partial, hasNA=hasNA, give.names=give.names)
 }
+frollvar = function(x, n, fill=NA, algo=c("fast","exact"), align=c("right","left","center"), na.rm=FALSE, has.nf=NA, adaptive=FALSE, partial=FALSE, give.names=FALSE, hasNA) {
+  froll(fun="var", x=x, n=n, fill=fill, algo=algo, align=align, na.rm=na.rm, has.nf=has.nf, adaptive=adaptive, partial=partial, hasNA=hasNA, give.names=give.names)
+}
+frollsd = function(x, n, fill=NA, algo=c("fast","exact"), align=c("right","left","center"), na.rm=FALSE, has.nf=NA, adaptive=FALSE, partial=FALSE, give.names=FALSE, hasNA) {
+  froll(fun="sd", x=x, n=n, fill=fill, algo=algo, align=align, na.rm=na.rm, has.nf=has.nf, adaptive=adaptive, partial=partial, hasNA=hasNA, give.names=give.names)
+}

inst/tests/froll.Rraw

@@ -1418,6 +1418,86 @@ test(6001.693, frollapply(FUN=median, adaptive=TRUE, c(1:2,NA), c(2,0,2)), c(NA,
 test(6001.694, frollapply(FUN=median, adaptive=TRUE, c(1:2,NA), c(2,0,2), na.rm=TRUE), c(NA,NA_integer_,2L))
 test(6001.695, frollapply(FUN=median, adaptive=TRUE, c(1:2,NA_real_), c(2,0,2), na.rm=TRUE, partial=TRUE), c(1,NA_real_,2))
 
+test(6001.711, frollvar(1:3, 0), c(NA_real_,NA_real_,NA_real_), options=c("datatable.verbose"=TRUE), output="window width of size 0")
+test(6001.712, frollvar(1:3, 0, fill=99), c(NA_real_,NA_real_,NA_real_))
+test(6001.713, frollvar(c(1:2,NA), 0), c(NA_real_,NA_real_,NA_real_))
+test(6001.714, frollvar(c(1:2,NA), 0, na.rm=TRUE), c(NA_real_,NA_real_,NA_real_))
+test(6001.715, frollvar(1:3, 0, algo="exact"), c(NA_real_,NA_real_,NA_real_), options=c("datatable.verbose"=TRUE), output="window width of size 0")
+test(6001.716, frollvar(c(1:2,NA), 0, algo="exact"), c(NA_real_,NA_real_,NA_real_))
+test(6001.717, frollvar(c(1:2,NA), 0, algo="exact", na.rm=TRUE), c(NA_real_,NA_real_,NA_real_))
+test(6001.718, frollvar(c(1:2,NA), 2), c(NA,0.5,NA), options=c("datatable.verbose"=TRUE), output="redirecting to frollvarExact")
+test(6001.721, frollvar(adaptive=TRUE, 1:3, c(2,0,2)), c(NA,NA,0.5), options=c("datatable.verbose"=TRUE), output="not implemented, fall back to")
+test(6001.722, frollvar(adaptive=TRUE, 1:3, c(2,0,2), fill=99), c(99,NA,0.5))
+test(6001.723, frollvar(adaptive=TRUE, c(1:2,NA), c(2,0,2)), c(NA_real_,NA_real_,NA_real_))
+test(6001.724, frollvar(adaptive=TRUE, c(1:2,NA), c(2,0,2), na.rm=TRUE), c(NA_real_,NA_real_,NA_real_))
+test(6001.7241, frollvar(adaptive=TRUE, c(1:2,NA), c(2,2,2), has.nf=FALSE), c(NA_real_,0.5,NA_real_), warning="used but non-finite values are present in input")
+test(6001.7242, frollvar(adaptive=TRUE, c(1:2,NA), c(2,2,2)), c(NA_real_,0.5,NA_real_), options=c("datatable.verbose"=TRUE), output="propagates NFs properply, no need to re-run")
+test(6001.7243, frollvar(adaptive=TRUE, c(1:2,NA), c(2,2,2), na.rm=TRUE), c(NA_real_,0.5,NA_real_), options=c("datatable.verbose"=TRUE), output="re-running with extra care for NFs")
+test(6001.725, frollvar(adaptive=TRUE, 1:3, c(2,0,2), algo="exact"), c(NA,NA,0.5))
+test(6001.726, frollvar(adaptive=TRUE, 1:3, c(2,0,2), fill=99, algo="exact"), c(99,NA,0.5))
+test(6001.727, frollvar(adaptive=TRUE, c(1:2,NA), c(2,0,2), algo="exact"), c(NA_real_,NA_real_,NA_real_))
+test(6001.728, frollvar(adaptive=TRUE, c(1:2,NA), c(2,0,2), algo="exact", na.rm=TRUE), c(NA_real_,NA_real_,NA_real_))
+test(6001.729, frollvar(adaptive=TRUE, c(1:2,NA), c(2,0,2), algo="exact", na.rm=TRUE, partial=TRUE), c(NA_real_,NA_real_,NA_real_))
+test(6001.730, frollvar(adaptive=TRUE, c(1:2,NA), c(2,0,2), fill=99, algo="exact", na.rm=TRUE), c(99,NA,NA))
+y = c(1e8+2.980232e-8, 1e8, 1e8, 1e8) # CLAMP0 test
+test(6001.731, frollvar(y, 3)[4L], 0)
+test(6001.732, frollsd(y, 3)[4L], 0)
+test(6001.733, frollvar(y, c(3,3,3,3), adaptive=TRUE)[4L], 0)
+test(6001.734, frollsd(y, c(3,3,3,3), adaptive=TRUE)[4L], 0)
+test(6001.781, frollapply(FUN=var, 1:3, 0), c(NA_real_,NA_real_,NA_real_))
+test(6001.782, frollapply(FUN=var, 1:3, 0, fill=99), c(NA_real_,NA_real_,NA_real_))
+test(6001.783, frollapply(FUN=var, c(1:2,NA), 0), c(NA_real_,NA_real_,NA_real_))
+test(6001.784, frollapply(FUN=var, c(1:2,NA), 0, na.rm=TRUE), c(NA_real_,NA_real_,NA_real_))
+test(6001.7910, frollapply(FUN=var, adaptive=TRUE, 1:3, c(2,0,2)), c(NA,NA,0.5))
+test(6001.7911, frollapply(FUN=var, adaptive=TRUE, list(1:3,2:4), c(2,0,2)), list(c(NA,NA,0.5), c(NA,NA,0.5)))
+test(6001.7912, frollapply(FUN=var, adaptive=TRUE, 1:3, list(c(2,0,2), c(0,2,0))), list(c(NA,NA,0.5), c(NA,0.5,NA)))
+test(6001.7913, frollapply(FUN=var, adaptive=TRUE, list(1:3,2:4), list(c(2,0,2), c(0,2,0))), list(c(NA,NA,0.5), c(NA,0.5,NA), c(NA,NA,0.5), c(NA,0.5,NA)))
+test(6001.792, frollapply(FUN=var, adaptive=TRUE, 1:3, c(2,0,2), fill=99), c(99,NA,0.5))
+test(6001.793, frollapply(FUN=var, adaptive=TRUE, c(1:2,NA), c(2,0,2)), c(NA_real_,NA_real_,NA_real_))
+test(6001.794, frollapply(FUN=var, adaptive=TRUE, c(1:2,NA), c(2,0,2), na.rm=TRUE), c(NA_real_,NA_real_,NA_real_))
+test(6001.795, frollapply(FUN=var, adaptive=TRUE, c(1:2,NA_real_), c(2,0,2), na.rm=TRUE, partial=TRUE), c(NA_real_,NA_real_,NA_real_))
+
+test(6001.810, frollsd(1:3, 0), c(NA_real_,NA_real_,NA_real_), options=c("datatable.verbose"=TRUE), output="frollsdFast: calling sqrt(frollvarFast(...))")
+test(6001.811, frollsd(1:3, 0), c(NA_real_,NA_real_,NA_real_), options=c("datatable.verbose"=TRUE), output="window width of size 0")
+test(6001.812, frollsd(1:3, 0, fill=99), c(NA_real_,NA_real_,NA_real_))
+test(6001.813, frollsd(c(1:2,NA), 0), c(NA_real_,NA_real_,NA_real_))
+test(6001.814, frollsd(c(1:2,NA), 0, na.rm=TRUE), c(NA_real_,NA_real_,NA_real_))
+test(6001.815, frollsd(1:3, 0, algo="exact"), c(NA_real_,NA_real_,NA_real_), options=c("datatable.verbose"=TRUE), output="window width of size 0")
+test(6001.816, frollsd(c(1:2,NA), 0, algo="exact"), c(NA_real_,NA_real_,NA_real_))
+test(6001.817, frollsd(c(1:2,NA), 0, algo="exact", na.rm=TRUE), c(NA_real_,NA_real_,NA_real_))
+test(6001.818, frollsd(c(1:2,NA), 2), c(NA,sqrt(0.5),NA), options=c("datatable.verbose"=TRUE), output="redirecting to frollvarExact")
+test(6001.8191, frollsd(1:3, 2, fill=99), c(99,sqrt(0.5),sqrt(0.5)))
+test(6001.8192, frollsd(1:3, 2, fill=99, algo="exact"), c(99,sqrt(0.5),sqrt(0.5)))
+test(6001.8193, frollsd(c(1:2,NA), 2, has.nf=FALSE), c(NA,sqrt(0.5),NA), warning="used but non-finite values are present in input")
+test(6001.8194, frollsd(c(NA,2:3), 2, has.nf=FALSE), c(NA,NA,sqrt(0.5)), warning="used but non-finite values are present in input")
+test(6001.8195, frollsd(c(NA,2:3), 2), c(NA,NA,sqrt(0.5)), options=c("datatable.verbose"=TRUE), output="skip non-finite inaware attempt and run with extra care")
+test(6001.8196, frollsd(c(NA,2:3), 2, has.nf=FALSE, algo="exact"), c(NA,NA,sqrt(0.5)), warning="used but non-finite values are present in input")
+test(6001.8197, frollsd(c(NA,2:3), 2, algo="exact", na.rm=TRUE), c(NA,NA,sqrt(0.5)), options=c("datatable.verbose"=TRUE), output="re-running with extra care for NF")
+test(6001.8201, frollsd(adaptive=TRUE, 1:3, c(2,2,2)), c(NA,sqrt(0.5),sqrt(0.5)), options=c("datatable.verbose"=TRUE), output="frolladaptivefun: algo 0 not implemented, fall back to 1")
+test(6001.8202, frollsd(adaptive=TRUE, 1:3, c(2,2,2)), c(NA,sqrt(0.5),sqrt(0.5)), options=c("datatable.verbose"=TRUE), output="frolladaptivesdExact: calling sqrt(frolladaptivevarExact(...))")
+test(6001.821, frollsd(adaptive=TRUE, 1:3, c(2,0,2)), c(NA,NA,sqrt(0.5)))
+test(6001.822, frollsd(adaptive=TRUE, 1:3, c(2,0,2), fill=99), c(99,NA,sqrt(0.5)))
+test(6001.823, frollsd(adaptive=TRUE, c(1:2,NA), c(2,0,2)), c(NA_real_,NA_real_,NA_real_))
+test(6001.824, frollsd(adaptive=TRUE, c(1:2,NA), c(2,0,2), na.rm=TRUE), c(NA_real_,NA_real_,NA_real_))
+test(6001.825, frollsd(adaptive=TRUE, 1:3, c(2,0,2), algo="exact"), c(NA,NA,sqrt(0.5)))
+test(6001.826, frollsd(adaptive=TRUE, 1:3, c(2,0,2), fill=99, algo="exact"), c(99,NA,sqrt(0.5)))
+test(6001.827, frollsd(adaptive=TRUE, c(1:2,NA), c(2,0,2), algo="exact"), c(NA_real_,NA_real_,NA_real_))
+test(6001.828, frollsd(adaptive=TRUE, c(1:2,NA), c(2,0,2), algo="exact", na.rm=TRUE), c(NA_real_,NA_real_,NA_real_))
+test(6001.829, frollsd(adaptive=TRUE, c(1:2,NA), c(2,0,2), algo="exact", na.rm=TRUE, partial=TRUE), c(NA_real_,NA_real_,NA_real_))
+test(6001.830, frollsd(adaptive=TRUE, c(1:2,NA), c(2,0,2), fill=99, algo="exact", na.rm=TRUE), c(99,NA,NA))
+test(6001.881, frollapply(FUN=sd, 1:3, 0), c(NA_real_,NA_real_,NA_real_))
+test(6001.882, frollapply(FUN=sd, 1:3, 0, fill=99), c(NA_real_,NA_real_,NA_real_))
+test(6001.883, frollapply(FUN=sd, c(1:2,NA), 0), c(NA_real_,NA_real_,NA_real_))
+test(6001.884, frollapply(FUN=sd, c(1:2,NA), 0, na.rm=TRUE), c(NA_real_,NA_real_,NA_real_))
+test(6001.8910, frollapply(FUN=sd, adaptive=TRUE, 1:3, c(2,0,2)), c(NA,NA,sqrt(0.5)))
+test(6001.8911, frollapply(FUN=sd, adaptive=TRUE, list(1:3,2:4), c(2,0,2)), list(c(NA,NA,sqrt(0.5)), c(NA,NA,sqrt(0.5))))
+test(6001.8912, frollapply(FUN=sd, adaptive=TRUE, 1:3, list(c(2,0,2), c(0,2,0))), list(c(NA,NA,sqrt(0.5)), c(NA,sqrt(0.5),NA)))
+test(6001.8913, frollapply(FUN=sd, adaptive=TRUE, list(1:3,2:4), list(c(2,0,2), c(0,2,0))), list(c(NA,NA,sqrt(0.5)), c(NA,sqrt(0.5),NA), c(NA,NA,sqrt(0.5)), c(NA,sqrt(0.5),NA)))
+test(6001.892, frollapply(FUN=sd, adaptive=TRUE, 1:3, c(2,0,2), fill=99), c(99,NA,sqrt(0.5)))
+test(6001.893, frollapply(FUN=sd, adaptive=TRUE, c(1:2,NA), c(2,0,2)), c(NA_real_,NA_real_,NA_real_))
+test(6001.894, frollapply(FUN=sd, adaptive=TRUE, c(1:2,NA), c(2,0,2), na.rm=TRUE), c(NA_real_,NA_real_,NA_real_))
+test(6001.895, frollapply(FUN=sd, adaptive=TRUE, c(1:2,NA_real_), c(2,0,2), na.rm=TRUE, partial=TRUE), c(NA_real_,NA_real_,NA_real_))
+
 # frollmedian
 rollmedian = function(x, k, na.rm=FALSE) {
   ans = rep(NA_real_, length(x))
@@ -2250,7 +2330,7 @@ rollfun = function(x, n, FUN, fill=NA_real_, na.rm=FALSE, nf.rm=FALSE, partial=F
   }
   ans
 }
-base_compare = function(x, n, funs=c("mean","sum","max","min","prod","median"), algos=c("fast","exact")) {
+base_compare = function(x, n, funs=c("mean","sum","max","min","prod","median","var","sd"), algos=c("fast","exact")) {
   num.step = 0.0001
   for (fun in funs) {
     for (na.rm in c(FALSE, TRUE)) {
@@ -2334,7 +2414,7 @@ base_compare(x, n)
 #### against zoo
 if (requireNamespace("zoo", quietly=TRUE)) {
   drollapply = function(...) as.double(zoo::rollapply(...)) # rollapply is not consistent in data type of answer, force to double
-  zoo_compare = function(x, n, funs=c("mean","sum","max","min","prod","median"), algos=c("fast","exact")) {
+  zoo_compare = function(x, n, funs=c("mean","sum","max","min","prod","median","var","sd"), algos=c("fast","exact")) {
     num.step = 0.0001
     #### fun, align, na.rm, fill, algo, partial
     for (fun in funs) {
@@ -2432,7 +2512,7 @@ arollfun = function(FUN, x, n, na.rm=FALSE, align=c("right","left"), fill=NA, nf
   }
   ans
 }
-afun_compare = function(x, n, funs=c("mean","sum","max","min","prod","median"), algos=c("fast","exact")) {
+afun_compare = function(x, n, funs=c("mean","sum","max","min","prod","median","var","sd"), algos=c("fast","exact")) {
   num.step = 0.0001
   #### fun, align, na.rm, fill, algo
   for (fun in funs) {

man/froll.Rd

@@ -10,13 +10,17 @@
 \alias{frollmin}
 \alias{frollprod}
 \alias{frollmedian}
+\alias{frollvar}
+\alias{frollsd}
 \alias{roll}
 \alias{rollmean}
 \alias{rollsum}
 \alias{rollmax}
 \alias{rollmin}
 \alias{rollprod}
 \alias{rollmedian}
+\alias{rollvar}
+\alias{rollsd}
 \title{Rolling functions}
 \description{
   Fast rolling functions to calculate aggregates on a sliding window. For a user-defined rolling function see \code{\link{frollapply}}. For "time-aware" (irregularly spaced time series) rolling function see \code{\link{frolladapt}}.
@@ -34,6 +38,10 @@
     na.rm=FALSE, has.nf=NA, adaptive=FALSE, partial=FALSE, give.names=FALSE, hasNA)
   frollmedian(x, n, fill=NA, algo=c("fast","exact"), align=c("right","left","center"),
     na.rm=FALSE, has.nf=NA, adaptive=FALSE, partial=FALSE, give.names=FALSE, hasNA)
+  frollvar(x, n, fill=NA, algo=c("fast","exact"), align=c("right","left","center"),
+    na.rm=FALSE, has.nf=NA, adaptive=FALSE, partial=FALSE, give.names=FALSE, hasNA)
+  frollsd(x, n, fill=NA, algo=c("fast","exact"), align=c("right","left","center"),
+    na.rm=FALSE, has.nf=NA, adaptive=FALSE, partial=FALSE, give.names=FALSE, hasNA)
 }
 \arguments{
   \item{x}{ Integer, numeric or logical vector, coerced to numeric, on which sliding window calculates an aggregate function. It supports vectorized input, then it needs to be a \code{data.table}, \code{data.frame} or a \code{list}, in which case a rolling function is applied to each column/vector. }
@@ -85,25 +93,26 @@
     \item \code{has.nf=TRUE} uses non-finite aware implementation straightaway.
     \item \code{has.nf=FALSE} uses faster implementation that does not support non-finite values. Then depending on the rolling function it will either:
     \itemize{
-      \item (\emph{mean, sum, prod}) detect non-finite, re-run non-finite aware.
+      \item (\emph{mean, sum, prod, var, sd}) detect non-finite, re-run non-finite aware.
       \item (\emph{max, min, median}) does not detect non-finites and may silently produce an incorrect answer.
     }
     In general \code{has.nf=FALSE && any(!is.finite(x))} should be considered as undefined behavior. Therefore \code{has.nf=FALSE} should be used with care.
   }
 }
 \section{Implementation}{
-  Each rolling function has 4 different implementations. First factor that decides which implementation is used is the \code{adaptive} argument (either \code{TRUE} or \code{FALSE}), see section below for details. Then for each of those two algorithms there are usually two implementations depending on the \code{algo} argument.
+  Most of the rolling functions have 4 different implementations. First factor that decides which implementation is used is the \code{adaptive} argument (either \code{TRUE} or \code{FALSE}), see section below for details. Then for each of those two algorithms there are usually two implementations depending on the \code{algo} argument.
   \itemize{
     \item \code{algo="fast"} uses \emph{"online"}, single pass, algorithm.
     \itemize{
       \item \emph{max} and \emph{min} rolling function will not do only a single pass but, on average, they will compute \code{length(x)/n} nested loops. The larger the window, the greater the advantage over the \emph{exact} algorithm, which computes \code{length(x)} nested loops. Note that \emph{exact} uses multiple CPUs so for a small window sizes and many CPUs it may actually be faster than \emph{fast}. However, in such cases the elapsed timings will likely be far below a single second.
       \item \emph{median} will use a novel algorithm described by \emph{Jukka Suomela} in his paper \emph{Median Filtering is Equivalent to Sorting (2014)}. See references section for the link. Implementation here is extended to support arbitrary length of input and an even window size. Despite extensive validation of results this function should be considered experimental. When missing values are detected it will fall back to slower \code{algo="exact"} implementation.
-      \item Not all functions have \emph{fast} implementation available. As of now adaptive \emph{max}, adaptive \emph{min} and adaptive \emph{median} do not have \emph{fast} implementation, therefore it will automatically fall back to \emph{exact} implementation. \code{datatable.verbose} option can be used to check that.
+      \item \emph{var} and \emph{sd} will use numerically stable \emph{Welford}'s online algorithm.
+      \item Not all functions have \emph{fast} implementation available. As of now, adaptive \emph{max}, \emph{min}, \emph{median}, \emph{var} and \emph{sd} do not have \emph{fast} adaptive implementation, therefore it will automatically fall back to \emph{exact} adaptive implementation. Similarly, non-adaptive fast implementations of \emph{median}, \emph{var} and \emph{sd} will fall back to \emph{exact} implementations if they detect any non-finite values in the input. \code{datatable.verbose} option can be used to check that.
     }
     \item \code{algo="exact"} will make the rolling functions use a more computationally-intensive algorithm. For each observation in the input vector it will compute a function on a rolling window from scratch (complexity \eqn{O(n^2)}).
     \itemize{
       \item Depeneding on the function, this algorithm may suffers less from floating point rounding error (the same consideration applies to base \code{\link[base]{mean}}).
-      \item In case of \emph{mean} (and possibly other functions in future), it will additionally make extra pass to perform floating point error correction. Error corrections might not be truly exact on some platforms (like Windows) when using multiple threads.
+      \item In case of \emph{mean}, it will additionally make an extra pass to perform floating point error correction. Error corrections might not be truly exact on some platforms (like Windows) when using multiple threads.
     }
   }
 }
@@ -158,6 +167,8 @@ frollmax(d, 3:4)
 frollmin(d, 3:4)
 frollprod(d, 3:4)
 frollmedian(d, 3:4)
+frollvar(d, 3:4)
+frollsd(d, 3:4)
 
 # partial=TRUE
 x = 1:6/2
@@ -179,6 +190,11 @@ frollsum(list(x=1:5, y=5:1), c(tiny=2, big=4), give.names=TRUE)
 frollmax(c(1,2,NA,4,5), 2)
 frollmax(c(1,2,NA,4,5), 2, has.nf=FALSE)
 
+# use verobse=TRUE for extra insight
+.op = options(datatable.verbose = TRUE)
+frollsd(c(1:5,NA,7:8), 4)
+options(.op)
+
 # performance vs exactness
 set.seed(108)
 x = sample(c(rnorm(1e3, 1e6, 5e5), 5e9, 5e-9))