Instantly share code, notes, and snippets.

# richarddmorey/logRepresentedReal.R Last active May 3, 2016

Representing real numbers using logarithms in R
 ########### # S4 class to represent real numbers as logarithms # Includes standard arithmetic operations +,-,*,/, and ^ # Richard D. Morey (richarddmorey@gmail.com) # November 2014 ########### # Compute log(exp(a) + exp(b)) logExpAplusExpB <- function(a, b) { a + log1plusExpA( b - a ) } # Compute log(1 + exp(a)) # From plogis.c log1plusExpA <- function(a) { if( a <= 18 ) return( log1p( exp(a) ) ) if( a > 33.3 ) return( a ) # else: 18.0 < x <= 33.3 : return( a + exp( -a ) ) } # Compute log(exp(a) - exp(b)) logExpAminusExpB <- function(a, b) { a + pexp( a - b, log.p=TRUE ) } setClass("logRepresentedReal", representation(modulo = "numeric" , sign ="integer") ) logRepresentedReal <- function( modulo, sign ){ new( "logRepresentedReal", modulo = modulo, sign = sign ) } setValidity("logRepresentedReal", function(object){ if( abs(object@sign) != 1 ) return("Sign must be 1 or -1.") return(TRUE) }) setMethod("show", signature("logRepresentedReal"), function(object){ # Use the function in the BayesFactor package to print # logarithmically-represented values nicely require(BayesFactor, quietly = TRUE) str = BayesFactor:::expString(object@modulo) sgn = "" if(object@sign < 0) sgn = "-" cat(paste(" ", sgn, str, " ", sep="")) }) setAs("logRepresentedReal", "numeric", function( from, to ){ as.double.logRepresentedReal(from) }) as.double.logRepresentedReal <- function(x){ x@sign * exp(x@modulo) } setMethod("abs", signature("logRepresentedReal"), function(x){ x@sign = 1L return(x) }) setMethod("sign", signature("logRepresentedReal"), function(x){ x@sign }) setMethod("-", signature("logRepresentedReal"), function(e1){ e1@sign = as.integer(-e1@sign) return(e1) }) setMethod('>', signature("logRepresentedReal", "logRepresentedReal"), function(e1, e2){ if( (e1@modulo == e2@modulo) & (e1@sign == e2@sign) ) return(FALSE) if(sign(e1) > sign(e2)) return(TRUE) if(sign(e1) < sign(e2)) return(FALSE) if(sign(e1) == 1) return(e1@modulo > e2@modulo) if(sign(e1) == -1) return(e1@modulo < e2@modulo) }) setMethod('<', signature("logRepresentedReal", "logRepresentedReal"), function(e1, e2){ if( (e1@modulo == e2@modulo) & (e1@sign == e2@sign) ) return(FALSE) return(!(e1>e2)) }) setGeneric("is.zero", function(x){ standardGeneric("is.zero") }) setMethod("is.zero", signature("logRepresentedReal"), function(x){ if( !is.infinite(x@modulo) ) return(FALSE) if(sign(x) != sign(x@modulo)) return(TRUE) return(FALSE) }) setMethod("is.infinite", signature("logRepresentedReal"), function(x){ if( !is.infinite(x@modulo) ) return(FALSE) if(sign(x) == sign(x@modulo)) return(TRUE) return(FALSE) }) setMethod("is.finite", signature("logRepresentedReal"), function(x){ return(!is.infinite(x)) }) setMethod('+', signature("logRepresentedReal", "logRepresentedReal"), function(e1, e2){ if(e1@sign == 1 & e2@sign == -1 ) return(e1 - -e2) if(e1@sign == -1 & e2@sign == -1 ) return( -(-e1 + -e2) ) if(e1@sign == -1 & e2@sign == 1 ) return(e2 - -e1) logsum = logExpAplusExpB(e1@modulo, e2@modulo) logRepresentedReal( logsum, 1L ) }) setMethod('+', signature("logRepresentedReal", "numeric"), function(e1, e2){ if( e2 == 0 ) return(e1) e2 = logRepresentedReal( log(abs(e2)), as.integer(sign(e2)) ) e1 + e2 }) setMethod('+', signature("numeric","logRepresentedReal"), function(e1, e2){ e2 + e1 }) setMethod('-', signature("logRepresentedReal", "logRepresentedReal"), function(e1, e2){ if(e1@sign == 1 & e2@sign == -1 ) return( e1 + -e2 ) if(e1@sign == -1 & e2@sign == -1 ) return( -abs(e1) + abs(e2) ) if(e1@sign == -1 & e2@sign == 1 ) return( -(e2 + -e1) ) if(e1 > e2){ logdiff = logExpAminusExpB(e1@modulo, e2@modulo) return(logRepresentedReal( logdiff, 1L )) }else if(e1 < e2){ logdiff = logExpAminusExpB(e2@modulo, e1@modulo) return(logRepresentedReal( logdiff, -1L )) }else{ return(logRepresentedReal( -Inf, 1L )) } }) setMethod('-', signature("logRepresentedReal", "numeric"), function(e1, e2){ if( e2 == 0 ) return(e1) e2 = logRepresentedReal( log(abs(e2)), as.integer(sign(e2)) ) e1 - e2 }) setMethod('-', signature("numeric","logRepresentedReal"), function(e1, e2){ -e2 + e1 }) setMethod('*', signature("numeric", "logRepresentedReal"), function(e1, e2){ sgn = sign(e1) * sign(e2) modu = log(abs(e1)) + e2@modulo logRepresentedReal( modu, as.integer(sgn) ) }) setMethod('*', signature("logRepresentedReal", "numeric"), function(e1, e2){ e2 * e1 }) setMethod('*', signature("logRepresentedReal", "logRepresentedReal"), function(e1, e2){ sgn = sign(e1) * sign(e2) modu = e1@modulo + e2@modulo logRepresentedReal( modu, as.integer(sgn) ) }) setMethod('/', signature("logRepresentedReal", "logRepresentedReal"), function(e1, e2){ sgn = sign(e1) * sign(e2) modu = e1@modulo - e2@modulo logRepresentedReal( modu, as.integer(sgn) ) }) setMethod('/', signature("numeric", "logRepresentedReal"), function(e1, e2){ if(e1 == 1){ e2@modulo = -e2@modulo return(e2) }else{ return(e1 * (1/e2)) } }) setMethod('/', signature("logRepresentedReal", "numeric"), function(e1, e2){ 1/(e2/e1) }) setMethod('^', signature("logRepresentedReal", "numeric"), function(e1, e2){ if( e2 == 0 ) return( logRepresentedReal( 0, 1L ) ) if( sign(e1)==1 | e2%%2 == 0 ){ sgn = 1L }else{ sgn = -1L } if(e2 < 0){ e2 = abs(e2) e1 = 1/e1 } modu = e1@modulo * e2 logRepresentedReal( modu, as.integer(sgn) ) })