From murdoch@dunc@n @ending from gm@il@com  Tue Sep 11 18:11:29 2018
From: murdoch@dunc@n @ending from gm@il@com (Duncan Murdoch)
Date: Tue, 11 Sep 2018 12:11:29 -0400
Subject: [Rd] Modification-proposal for %% (modulo) when supplied with
 double
In-Reply-To: <1BC486D5-5B53-4C23-B422-877B00EF91B5@dans.knaw.nl>
References: <1BC486D5-5B53-4C23-B422-877B00EF91B5@dans.knaw.nl>
Message-ID: <15346010-f2c5-8aa9-7b93-5e24c0692c32@gmail.com>

On 11/09/2018 11:23 AM, Emil Bode wrote:
> Hi all,
> 
> 
> 
> Could we modify the "%%" (modulo)-operator to include some tolerance for rounding-errors when supplied with doubles?
> 
> It's not much work (patch supplied on the bottom), and I don't think it would break anything, only if you were really interested in analysing rounding differences.
> 
> Any ideas about implementing this and overwriting base::`%%`, or would we want another method (as I've done for the moment)?

I think this is a bad idea.  Your comments say "The 
\code{\link[base:Arithmetic]{`\%\%`}} operator calculates the modulo, 
but sometimes has rounding errors, e.g. "\code{(9.1/.1) \%\% 1}" gives ~ 
1, instead of 0."

This is false.  The %% calculation is exactly correct.  The rounding 
error happened in your input:  9.1/0.1 is not equal to 91, it is a 
little bit less:

 > options(digits=20)
 > 9.1/.1
[1] 90.999999999999985789

And %% did not return 1, it returned the correct value:

 > (9.1/.1) %% 1
[1] 0.99999999999998578915

So it makes no sense to change %%.

You might argue that the division 9.1/.1 is giving the wrong answer, but 
in fact that answer is correct too.  The real problem is that in double 
precision floating point the numbers 9.1 and .1 can't be represented 
exactly.  This is well known, it's in the FAQ (question 7.31).

Duncan Murdoch

> 
> 
> 
> Background
> 
> I was writing some code where something has to happen at a certain interval, with progress indicated, something like this:
> 
> interval <- .001
> 
> progress <- .1
> 
> for(i in 1:1000*interval) {myFun(i); Sys.sleep(interval); if(i %% progress, 0))) cat(i, '\n')}
> 
> without interval and progress being known in advance. I could work around it and make i integer, or do something like
> 
> isTRUE(all.equal(i %% progress,0)) || isTRUE(all.equal(i %% progress, progress),
> 
> but I think my code is clearer as it is. And I like the idea behind all.equal: we want double to approximately identical.
> 
> 
> 
> So my patch (with roxygen2-markup):
> 
> #' Modulo-operator with near-equality
> 
> #'
> 
> #' The \code{\link[base:Arithmetic]{`\%\%`}} operator calculates the modulo, but sometimes has rounding errors, e.g. "\code{(9.1/.1) \%\% 1}" gives ~ 1, instead of 0.\cr
> 
> #' Comparable to what all.equal does, this operator has some tolerance for small rounding errors.\cr
> 
> #' If the answer would be equal to the divisor within a small tolerance, 0 is returned instead.
> 
> #'
> 
> #' For integer x and y, the normal \%\%-operator is used
> 
> #'
> 
> #' @usage `\%mod\%`(x, y, tolerance = sqrt(.Machine$double.eps))
> 
> #' x \%mod\% y
> 
> #' @param x,y numeric vectors, similar to those passed on to \%\%
> 
> #' @param tolerance numeric, maximum difference, see \code{\link[base]{all.equal}}. The default is ~ \code{1.5e-8}
> 
> #' @return identical to the result for \%\%, unless the answer would be really close to y, in which case 0 is returned
> 
> #' @note To specify tolerance, use the call \code{`\%mod\%`(x,y,tolerance)}
> 
> #' @note The precedence for \code{\%mod\%} is the same as that for \code{\%\%}
> 
> #'
> 
> #' @name mod
> 
> #' @rdname mod
> 
> #'
> 
> #' @export
> 
> `%mod%` <- function(x,y, tolerance = sqrt(.Machine$double.eps)) {
> 
>    stopifnot(is.numeric(x), is.numeric(y), is.numeric(tolerance),
> 
>              !is.na(tolerance), length(tolerance)==1, tolerance>=0)
> 
>    if(is.integer(x) && is.integer(y)) {
> 
>      return(x %% y)
> 
>    } else {
> 
>      ans <- x %% y
> 
>      return(ifelse(abs(ans-y)<tolerance | abs(ans)<tolerance, 0, ans))
> 
>    }
> 
> }
> 
> 
> 
> Best regards,
> 
> Emil Bode
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>



[Rd] Modification-proposal for %% (modulo) when supplied with double
Duncan Murdoch