\documentclass[nojss]{jss} \usepackage{dsfont} \usepackage{bbm} \usepackage{amsfonts} \usepackage{wasysym} \usepackage{wrapfig} \author{Robin K. S. Hankin\\Auckland University of Technology} \title{A step-by-step guide to writing a simple package that uses \proglang{S4} methods: a ``hello world'' example} %\VignetteIndexEntry{Brobdingnag: a ``hello world'' package using S4 methods} \Plainauthor{Robin K. S. Hankin} \Plaintitle{Brobdingnagian numbers in S4} \Shorttitle{Brobdingnagian numbers in S4} \Abstract{ This vignette shows how to use \proglang{S4} methods to create a simple package. The other vignette shows how to use the package for solving problems involving very large numbers; it is based on~\cite{Rnews:Hankin:2007}. } \Keywords{\proglang{S4} methods, Brobdingnag, \proglang{R}} \Plainkeywords{S4 methods, Brobdingnag, R} \Address{ Robin K. S. Hankin\\ Auckland University of Technology\\ E-mail: \email{hankin.robin@gmail.com} } %% need no \usepackage{Sweave.sty} \SweaveOpts{echo=TRUE} \begin{document} \newsymbol\leqslant 1336 \section{Introduction} \setlength{\intextsep}{0pt} \begin{wrapfigure}{r}{0.2\textwidth} \begin{center} \includegraphics[width=1in]{\Sexpr{system.file("help/figures/Brobdingnag.png",package="Brobdingnag")}} \end{center} \end{wrapfigure} This vignette proves that it is possible for a `normal' person\footnote{That is, someone without super powers (such as might manifest themselves after being bitten by a radioactive member of \proglang{R}-core, for example)} to write a package using \proglang{S4} methods. It gives a step-by-step guide to creating a package that contains two \proglang{S4} classes (brobs and glubs) and a bunch of basic utilities for manipulating them. This document focuses on the \proglang{S4} aspects of the package. For an overview of the mathematical properties of Brobdingnagian numbers, their potential and their limitations, see the {\tt .Rd} files and~\cite{Rnews:Hankin:2007}. If you like this vignette and package, and find it useful, let me know. If there is anything wrong with it, let me know. My current thinking is that \proglang{S4} methods are more flexible and scalablee than \proglang{S3}, but much more work to set up (many of my packages use \proglang{S3} methods which I found to be perfectly adequate for my needs). Reasons for using \proglang{S4} might include a package having a large number of object classes that have a complicated hierarchical structure, or a complicated set of methods that interact with the object classes in a complicated manner. In the package, brobs are dealt with in {\tt brob.R}, and glubs are treated in {\tt glub.R} which appropriately generalizes all the {\tt brob} functionality. This document could not have been prepared (and should not be read) without consulting the following resources: \begin{itemize} \item John M. Chambers (1998), {\it Programming with Data}. New York: Springer, ISBN 0-387-98503-4 (The Green Book). \item W. N. Venables and B. D. Ripley (2000), {\it S Programming}. Springer, ISBN 0-387-98966-8. \item John Chambers (2006). {\tt How \proglang{S4} methods work} (available on CRAN). \end{itemize} \subsection{Overview} The idea of the \pkg{Brobdingnag} package is simple: the IEEE representation for floating point numbers cannot represent numbers larger than about~$1.8\times 10^{308}$. The package represents a number by the natural logarithm of its magnitude, and also stores a Boolean value indicating its sign. Objects so stored have class {\tt brob}; complex numbers may be similarly represented and have class {\tt glub}. With this scheme, multiplication is easy but addition is hard. The basic identity is: \[ \log(e^x+e^y) = \left\{ \begin{array}{cc} x+\log\left(1+e^{y-x}\right)\qquad &\mbox{if~$x>y$}\\ y+\log\left(1+e^{x-y}\right)\qquad &\mbox{otherwise} \end{array} \right. \] In practice this gets more complicated as one has to keep track of the sign; and special dispensation is needed for zero ($=e^{-\infty}$). One can thus deal with numbers up to about~$e^{1.9\times 10^{308}}\simeq 10^{7.8\times 10^{307}}$, although at this outer limit accuracy is pretty poor. \section{Class definition} The first thing we need to do is to define the {\tt brob} class. This uses the {\tt setClass()} function: <>= <>= setClass("swift", representation = "VIRTUAL" ) setClass("brob", representation = representation(x="numeric",positive="logical"), prototype = list(x=numeric(),positive=logical()), contains = "swift" ) @ It is simpler to ignore the first call to {\tt setClass()} here; for reasons that will become apparent when discussing the {\tt c()} function, one needs a virtual class that contains class brob and glub. To understand virtual classes, see section~\ref{logicSection}. The second call to {\tt setClass()} is more germane. Let's take this apart, argument by argument. The first argument, {\tt representation}, specifies ``the slots that the new class should have and/or other classes that this class extends. Usually a call to the `representation' function''. The helppage for {\tt representation} gives a few further details. Thus this argument specifies two `slots': one for the value and one for the sign. These are specified to be numeric and logical (NB: not Boolean) respectively. The second argument, {\tt prototype}, specifies default data for the slots. This kicks in when defining a zero-length brob; an example would be extracting {\tt x[FALSE]} where {\tt x} is a brob. The third argument, {\tt contains}, tells \proglang{R} that class {\tt swift} (which was specified to be virtual), has {\tt brob} as a subclass. We will need this later when we start to deal with {\tt glub}s, which are also a subclass of {\tt swift}. Let's use it: <>= new("brob",x=1:10,positive=rep(TRUE,10)) @ Notes: \begin{itemize} \item Function {\tt new()} is the {\em only} way to create objects of class brob. So, any object of class brob {\em must} have been created with function {\tt new()}. This is part of what the ``formal'' tag for \proglang{S4} means\footnote{Compare \proglang{S3}, in which I can say {\tt a <- 1:10; class(a) <- "lilliput"}}. \item Function {\tt new()} requires its arguments to be named, and no partial argument matching is performed. \item Function {\tt new()} is not intended for the user. It's too picky and difficult. To create new brobs, we need some more friendly functions---{\tt as.brob()} and {\tt brob()}---discussed below. \item There is, as yet, no print method, so the form of the object printed to the screen is less than ideal. \end{itemize} \subsection{Validity methods} Now, an optional step is to define a function that tests whether the arguments passed to {\tt new()} are acceptable. As it stands, the following code: <>= new("brob",x=1:10,positive=c(TRUE,FALSE,FALSE)) @ \noindent will not return an error, but is not acceptable because the arguments are different lengths (and will not recycle neatly). So, we define a validity method: <>= .Brob.valid <- function(object){ len <- length(object@positive) if(len != length(object@x)){ return("length mismatch") } else { return(TRUE) } } @ Advice on the form of the validity-testing function---here {\tt .Brob.valid()}---is given in the help page for {\tt setValidity}: ``The method should be a function of one object that returns `TRUE' or a description of the non-validity''. Examples are given in section 7.1.6 of the Green Book. In this package, I define a whole bunch of functions whose name starts with {\tt .Brob.}; these are internal and not intended for the user. They are also not documented. So now we have a function, {\tt .Brob.valid()}, that checks whether its argument has slots of the same length. We need to tell \proglang{R} that this function should be invoked every time a {\tt brob} is created. Function {\tt setValidity()} does this: <>= setValidity("brob", .Brob.valid) @ Thus, from now on [ie after the above call to {\tt setValidity()}], when calling {\tt new("brob", ...)} the two arguments {\tt x} and {\tt positive} must be the same length: recycling is not carried out. Functions like {\tt .Brob.valid()} that are user-unfriendly all have names beginning with {\tt .Brob}. These functions are there to help the organization of the package and are not intended to be used by the end-user. Clever, user-friendly operations such as recycling are carried out in the more user-friendly functions such as {\tt as.brob()}. If one were to call {\tt new()} with arguments of differing lengths, as in\\ \\ {\tt new("brob",x=1:10,positive=TRUE) } \\ \\ then {\tt new()} would report the error message in function {\tt .Brob.valid()}, because the {\tt positive} argument had length~1 and the {\tt x} was length~10; and the validity method {\tt .Brob.valid()} requires both arguments to be the same length\footnote{Placing the above call in {\tt try()} and showing the error explicitly would cause the package to fail {\tt R CMD check}.}. So now {\tt new()} works, but isn't exactly user-friendly: often one would want the above call to recycle the second argument to length 10 to match the first. This deficiency is remedied in the next section. \section{Basic user-friendly functions to create brobs} The basic, semi user-friendly function for creating brobs is {\tt brob()}: <>= "brob" <- function(x=double(),positive){ if(missing(positive)){ positive <- rep(TRUE,length(x)) } if(length(positive)==1){ positive <- rep(positive,length(x)) } new("brob",x=as.numeric(x),positive=positive) } @ Thus {\tt brob(x)} will return a number formally equal to~$e^x$. Function {\tt brob()} does helpful things like assuming the user desires positive numbers; it also carries out recycling: <>= brob(1:10,FALSE) @ Note that {\tt brob()} isn't exactly terribly user-friendly: it's confusing. {\tt brob(5)} returns a number formally equal to~$e^5$, not~$5$. This is documented in the help page, where the user is encouraged to use function {\tt as.brob()} instead. \section[Testing for brobs: an ``is.brob()'' function]{Testing for brobs: an {\tt is.brob()} function} Function {\tt is()} will test for an object being a brob: <>= is(brob(1:5),"brob") @ (see {\tt help(is)} for more details) but a small package like this, with only brobs and glubs to consider, could benefit from an \proglang{S3}-style function {\tt is.brob()}. This is easy to define: <>= is.brob <- function(x){is(x,"brob")} is.glub <- function(x){is(x,"glub")} @ now the user can just type {\tt is.brob(x)} to find out whether an object is a brob\footnote{This approach is fine for a tiny package like Brobdingnag, with only two or three classes. However, in the context of a more complicated package such as {\tt Matrix}, which uses dozens of different classes in a complicated hierarchical structure, one might prefer to type {\tt is(x,"dpoMatrix-class")} rather than define a plethora of functions along the lines of {\tt is.dpoMatrix-class()}.}. We also define an {\tt is.glub()} function similarly. So now we can check for objects being brobs and glubs. \section[Coercion: an ``as.brob()'' function]{Coercion: an {\tt as.brob()} function} Next, some ways to coerce objects to class brob: <>= "as.brob" <- function(x){ if(is.brob(x)){ return(x) } else if(is.complex(x)) { warning("imaginary parts discarded") return(Recall(Re(x))) } else if(is.glub(x)){ warning("imaginary parts discarded") return(Re(x)) } else { return(brob(log(abs(x)), x>=0)) } } @ So now we can coerce objects of various classes to brobs\footnote{The recommended way, appropriate for a complicated package such as Matrix, would be to execute {\tt setAs("numeric", "brob", .Brob.numeric\_to\_brob)} and {\tt setAs("complex", "brob", .Brob.complex\_to\_brob)}. Then, if {\tt x} is numeric, {\tt as(x,"brob")} would return the appropriate brob via {\tt .Brob.numeric\_to\_brob()}; we could then make {\tt as.brob()} a generic with {\tt setGeneric()} and define methods for it. Doing it this way would save function {\tt as.brob()} having to test for its argument being a brob [carried out in the first few lines of {\tt as.brob()}] because {\tt as()} has such a test built in, implicitly.}. The {\em only} way to create a brob is to use {\tt new()}, and the {\em only} function that calls this is {\tt brob()}. And {\tt as.brob()} calls this. Note the user-friendliness of {\tt as.brob()}. It takes numerics, brobs, and glubs (which give a warning). Check {\tt as.brob()}: <>= as.brob(1:10) @ \section[Coercion: an ``as.numeric()'' function]{Coercion: an {\tt as.numeric()} function} Now we need some methods to coerce brobs to numeric. This is a two-stage process. <>= <>= setAs("brob", "numeric", function(from){ out <- exp(from@x) out[!from@positive] <- -out[!from@positive] return(out) } ) @ This call to {\tt setAs()} makes {\tt as(x,"numeric")} carry out the function passed as the third argument when given a brob. But in this package, the user isn't supposed to type {\tt as(x,"numeric")}: the user is supposed to type {\tt as.numeric(x)}\footnote{Users can be expected to be familiar with functions such as {\tt as.numeric()} and {\tt as.complex()}, which is why the Brobdingnag package recommends this form. However, this might not be appropriate for a more complicated package such as Matrix because, to quote Martin Maechler, ``it seems very ugly to define more than very few {\tt as.FOO()} methods, and very natural to work with {\tt as(*, "FOO")} [constructions]''. Of course, there's nothing to stop a user typing {\tt as(x,"numeric")} if they wish.}. To accomplish this, we have to tell \proglang{R} that function {\tt as.numeric()} should execute {\tt as(x,"numeric")} when given a {\tt brob}. This is done by calling function {\tt setMethod()}: <>= setMethod("as.numeric",signature(x="brob"),function(x){as(x,"numeric")}) @ We similarly need to make {\tt as.complex()} work for brobs: <>= setAs("brob", "complex", function(from){ return(as.numeric(from)+ 0i) } ) setMethod("as.complex",signature(x="brob"),function(x){as(x,"complex")}) @ We'll need similar methods for glubs too. Better check: <>= x <- as.brob(1:4) x as.numeric(x) @ So that works. \section{Print methods} Print methods are not strictly necessary, but make using the package much easier. First, a helper function: <>= .Brob.print <- function(x, digits=5){ noquote( paste(c("-","+")[1+x@positive],"exp(",signif(x@x,digits),")",sep="")) } @ Then an \proglang{S3} method: <>= print.brob <- function(x, ...){ jj <- .Brob.print(x, ...) print(jj) return(invisible(jj)) } @ And finally a call to {\tt setMethod()}: <>= setMethod("show", "brob", function(object){print.brob(object)}) @ This two-stage methodology is recommended in the Venables and Ripley. The {\tt .Brob.print()} function does the hard work. Example of it in use: <>= as.brob(1:4) @ See how the brob object is printed out nicely, and with no special effort required of the user. \section{Get and Set methods} To be anal retentive about things, one should define {\tt C++} style accessor functions as follows: <>= <>= setGeneric("getX",function(x){standardGeneric("getX")}) setGeneric("getP",function(x){standardGeneric("getP")}) setMethod("getX","brob",function(x){x@x}) setMethod("getP","brob",function(x){x@positive}) @ but in practice I just use {\tt @} to access the slots. These are just here for good form's sake. \section{Length} Now a length: <>= <>= setMethod("length","brob",function(x){length(x@x)}) @ \section{Extracting elements of a vector} Next thing is to define some methods for extraction. This is done with {\tt setMethod()} for extraction, and {\tt setReplaceMethod()} for replacement. <>= <>= setMethod("[", "brob", function(x, i, j, drop){ if(!missing(j)){ warning("second argument to extractor function ignored") } brob(x@x[i], x@positive[i]) } ) @ See how the third argument to {\tt setMethod()} is a function whose arguments are the same as those to {\tt "["() %] }. Argument {\tt j} {\em must} be there otherwise one gets a signature error. I've put in a warning if a second argument that might be interpreted as {\tt j} is given. Now a method for replacement. This is a call to {\tt setReplaceMethod()}: <>= <>= setReplaceMethod("[",signature(x="brob"), function(x,i,j,value){ if(!missing(j)){ warning("second argument to extractor function ignored") } jj.x <- x@x jj.pos <- x@positive if(is.brob(value)){ jj.x[i] <- value@x jj.pos[i] <- value@positive return(brob(x=jj.x,positive=jj.pos)) } else { x[i] <- as.brob(value) return(x) } } ) @ See how the replacement function tests for the replacement value being a brob and acts accordingly. \section[Concatenation function ``cbrob()'']{Concatenation function {\tt cbrob()}} \label{cbrob} It is not possible to make {\tt c()} behave as expected for brobs\footnote{The ideas in this section are entirely due to John Chambers, who kindly replied to a question of mine on the \proglang{R}-devel email list} (that is, if any of its arguments are brobs, to coerce all its arguments to brobs and then concatenate). However, it is possible to define a function {\tt cbrob()} that does the job. This has to be done in several stages. First we define another user-unfriendly helper function {\tt .Brob.cPair()} which takes two arguments, coerces them to brobs, and concatenates them: <<.Brob.cPair>>= .Brob.cPair <- function(x,y){ x <- as.brob(x) y <- as.brob(y) brob(c(x@x,y@x),c(x@positive,y@positive)) } @ This is just {\tt c()} for the two slots separately. The idea is that function {\tt .Brob.cPair()} takes two arguments; both are coerced to brobs and it returns the concatenated vector of brobs. Now, we need to set up a (user-unfriendly) generic function {\tt .cPair()}: <>= <>= setGeneric(".cPair", function(x,y){standardGeneric(".cPair")}) @ Function {\tt .cPair()} is not substantive (sic): it exists purely in order to be a generic function that dispatches to {\tt .Brob.cPair()}. Now we use {\tt setMethod()} to organize the dispatch: <>= <>= setMethod(".cPair", c("brob", "brob"), function(x,y){.Brob.cPair(x,y)}) setMethod(".cPair", c("brob", "ANY"), function(x,y){.Brob.cPair(x,as.brob(y))}) setMethod(".cPair", c("ANY", "brob"), function(x,y){.Brob.cPair(as.brob(x),y)}) setMethod(".cPair", c("ANY", "ANY"), function(x,y){c(x,y)}) @ The four calls are necessary for the four different signatures that might be encountered. Note the {\tt ANY} class in the second, third, and fourth call. Thus if someone wants to write a new class of object (a lugg, say), and wants to concatenate luggs with a brob, this will work provided that they use {\tt setAs()} to make {\tt as.brob()} coerce correctly for lugg objects. The method used here allows this to be done without any changes to the Brobdingnag package. The final stage is the definition of {\tt cbrob()}, a user-friendly wrapper for the above stuff: <>= "cbrob" <- function(x, ...) { if(nargs()<3) .cPair(x,...) else .cPair(x, Recall(...)) } @ Note the recursive definition. If {\tt cbrob()} is called with {\em any} set of arguments that include a brob anywhere, this will result in the whole lot being coerced to brobs [by {\tt .Brob.cPair()}]. Which is what we want (although glubs will require more work). Just test this: <>= a <- 1:3 b <- as.brob(1e100) cbrob(a,a,b,a) @ So it worked: everything was coerced to a brob because of the single object of class brob in the call. \section{Maths} The math group generic functions are set with function {\tt setMethod()}. But, before this can be called, function {\tt sqrt()} needs a specific brob method: <>= <>= setMethod("sqrt","brob", function(x){ brob(ifelse(x@positive,x@x/2, NaN),TRUE) } ) @ Just check that: <>= sqrt(brob(4)) @ With these out of the way we can use {\tt setMethod()} to define the appropriate functions in the math group generic: <>= <>= setMethod("Math", "brob", function(x){ switch(.Generic, abs = brob(x@x), log = { out <- x@x out[!x@positive] <- NaN out }, exp = brob(x), cosh = {(brob(x) + brob(-x))/2}, sinh = {(brob(x) - brob(-x))/2}, acos =, acosh =, asin =, asinh =, atan =, atanh =, cos =, sin =, tan =, tanh =, trunc = callGeneric(as.numeric(x)), lgamma =, cumsum =, gamma =, ceiling=, floor = as.brob(callGeneric(as.numeric(x))), stop(paste(.Generic, "not allowed on Brobdingnagian numbers")) ) } ) @ See how the third argument to {\tt setMethod()} is a function. This function has access to {\tt .Generic}, in addition to {\tt x} and uses it to decide which operation to perform. See how functions {\tt acos()} to {\tt trunc()} just drop through to {\tt callGeneric(as.numeric(x))}. See also the method for {\tt log()}, which uses facts about brobs not known to \proglang{S4}. Just a quick check: <>= sin(brob(4)) @ So that works. \section{Operations} Now we need to make sure that {\tt brob(1) + brob(3)} works: the operations {\tt +, -, *, /} must work as expected. This is hard. First step: define some user-unfriendly functions that carry out the operations. For example, function {\tt .Brob.negative()} simply returns the negative of a brob. These functions are not for the user. <<.brob.arithstuff>>= .Brob.negative <- function(e1){ brob(e1@x,!e1@positive) } .Brob.ds <- function(e1,e2){ xor(e1@positive,e2@positive) } .Brob.add <- function(e1,e2){ e1 <- as.brob(e1) e2 <- as.brob(e2) jj <- rbind(e1@x,e2@x) x1 <- jj[1,] x2 <- jj[2,] out.x <- double(length(x1)) jj <- rbind(e1@positive,e2@positive) p1 <- jj[1,] p2 <- jj[2,] out.pos <- p1 ds <- .Brob.ds(e1,e2) ss <- !ds out.x[ss] <- pmax(x1[ss],x2[ss]) + log1p(+exp(-abs(x1[ss]-x2[ss]))) out.x[ds] <- pmax(x1[ds],x2[ds]) + log1p(-exp(-abs(x1[ds]-x2[ds]))) out.x[ (x1 == -Inf) & (x2 == -Inf)] <- -Inf out.pos <- p1 out.pos[ds] <- xor((x1[ds] > x2[ds]) , (!p1[ds]) ) return(brob(out.x,out.pos)) } .Brob.mult <- function(e1,e2){ e1 <- as.brob(e1) e2 <- as.brob(e2) return(brob(e1@x + e2@x, !.Brob.ds(e1,e2))) } .Brob.power <- function(e1,e2){ stopifnot(is.brob(e1) | is.brob(e2)) if(is.brob(e2)){ return(brob(log(e1) * brob(e2@x), TRUE)) } else { s <- as.integer(2*e1@positive-1) return(brob(e1@x*as.brob(e2), (s^as.numeric(e2))>0)) } } .Brob.inverse <- function(b){brob(-b@x,b@positive)} @ Note the complexity of {\tt .Brob.add()}. This is hard because logs are good at multiplying but bad at adding [{\tt ss} and {\tt ds} mean ``same sign'' and ``different sign'' respectively]. The first step is to make sure that the unary operators {\tt +} and {\tt -} work. We do this by a call to {\tt setMethod()}: <>= <>= setMethod("Arith",signature(e1 = "brob", e2="missing"), function(e1,e2){ switch(.Generic, "+" = e1, "-" = .Brob.negative(e1), stop(paste("Unary operator", .Generic, "not allowed on Brobdingnagian numbers")) ) } ) @ Note the second argument is {\tt signature(e1 = "brob", e2="missing")}: this effectively restricts the scope to unary operators. The {\tt switch()} statement only allows the + and the -.\ Check: <>= -brob(5) @ So that works. Next step, another user-unfriendly helper function that does the dirty work: <>= .Brob.arith <- function(e1,e2){ switch(.Generic, "+" = .Brob.add (e1, e2), "-" = .Brob.add (e1, .Brob.negative(as.brob(e2))), "*" = .Brob.mult (e1, e2), "/" = .Brob.mult (e1, .Brob.inverse(as.brob(e2))), "^" = .Brob.power(e1, e2), stop(paste("binary operator \"", .Generic, "\" not defined for Brobdingnagian numbers")) ) } @ And now we can call {\tt setMethod()}: <>= setMethod("Arith", signature(e1 = "brob", e2="ANY"), .Brob.arith) setMethod("Arith", signature(e1 = "ANY", e2="brob"), .Brob.arith) setMethod("Arith", signature(e1 = "brob", e2="brob"), .Brob.arith) @ Better check it: <>= 1e100 + as.brob(10)^100 @ \section{Comparison} This is pretty much the same as the others. First, some user-unfriendly helper functions: <>= .Brob.equal <- function(e1,e2){ (e1@x==e2@x) & (e1@positive==e2@positive) } .Brob.greater <- function(e1,e2){ jj.x <- rbind(e1@x,e2@x) jj.p <- rbind(e1@positive,e2@positive) ds <- .Brob.ds(e1,e2) ss <- !ds greater <- logical(length(ss)) greater[ds] <- jj.p[1,ds] greater[ss] <- jj.p[1,ss] & (jj.x[1,ss] > jj.x[2,ss]) return(greater) } @ These are the fundamental ones. We can now define another all-encompassing user-unfriendly function: <>= ".Brob.compare" <- function(e1,e2){ e1 <- as.brob(e1) e2 <- as.brob(e2) switch(.Generic, "==" = .Brob.equal(e1,e2), "!=" = !.Brob.equal(e1,e2), ">" = .Brob.greater(e1,e2), "<" = !.Brob.greater(e1,e2) & !.Brob.equal(e1,e2), ">=" = .Brob.greater(e1,e2) | .Brob.equal(e1,e2), "<=" = !.Brob.greater(e1,e2) | .Brob.equal(e1,e2), stop(paste(.Generic, "not supported for Brobdingnagian numbers")) ) } @ See how this function coerces both arguments to brobs. Now the call to {\tt setMethod()}: <>= <>= setMethod("Compare", signature(e1="brob", e2="ANY" ), .Brob.compare) setMethod("Compare", signature(e1="ANY" , e2="brob"), .Brob.compare) setMethod("Compare", signature(e1="brob", e2="brob"), .Brob.compare) @ Better check: <>= as.brob(10) < as.brob(11) as.brob(10) <= as.brob(10) @ So that works. \section{Logic} \label{logicSection} (The material in this section works in \proglang{R}-2.4.1, but not \proglang{R}-2.4.0). First a helper function: <>= .Brob.logic <- function(e1,e2){ stop("No logic currently implemented for Brobdingnagian numbers") } @ Now the calls to {\tt setMethod()}: <>= <>= setMethod("Logic",signature(e1="swift",e2="ANY"), .Brob.logic) setMethod("Logic",signature(e1="ANY",e2="swift"), .Brob.logic) setMethod("Logic",signature(e1="swift",e2="swift"), .Brob.logic) @ Note that the signatures specify {\tt swift} objects, so that glubs will be handled correctly too in one fell swoop. Note the third call to {\tt setMethod()}: without this, a call to {\tt Logic()} with signature {\tt c("swift","swift")} would be ambiguous as it might be interpreted as {\tt c("swift","ANY")} or {\tt c("ANY","swift")}. Here, class {\tt swift} extends {\tt brob} and {\tt glub}; so both brobs and glubs {\em are} {\tt swift} objects. But no object is a ``pure'' {\tt swift}; it's either a brob or a glub. This is useful here because I might dream up some new class of objects that are ``like'' brobs in some way (for example, a class of objects whose elements are quaternions with Brobdingnagian components) and it would be nice to specify behaviour that is generic to brobs and glubs and the new class of objects too. \section{Miscellaneous generics} We now have to tell \proglang{R} that certain functions are to be considered generic. The functions are {\tt max(), min(), range(), prod()}, and {\tt sum()}. The help page for (eg) {\tt max()} specifies that the arguments must be numeric, and brobs aren't numeric. In versions of \proglang{R} prior to 2.6-0 (I think), {\tt log()} needs to be made a generic with {\tt setGeneric()}. But, in versions 2.6-0 and above, all the group generics (including {\tt log}) are primitive, which means that the generic function is implicit and cannot be changed. This applies to the other group generics too ({\tt max}, {\tt min}, {\tt prod}, {\tt range} {\tt sum}). So to work with both types of \proglang{R}, one needs to check whether or not {\tt log} (or {\tt max} or whatever) is generic before calling {\tt setGeneric()}. <>= if(!isGeneric("log")){ setGeneric("log",group="Math") } @ <>= <>= if(!isGeneric("sum")){ setGeneric("max", function(x, ..., na.rm = FALSE) { standardGeneric("max") }, useAsDefault = function(x, ..., na.rm = FALSE) { base::max(x, ..., na.rm = na.rm) }, group = "Summary") setGeneric("min", function(x, ..., na.rm = FALSE) { standardGeneric("min") }, useAsDefault = function(x, ..., na.rm = FALSE) { base::min(x, ..., na.rm = na.rm) }, group = "Summary") setGeneric("range", function(x, ..., na.rm = FALSE) { standardGeneric("range") }, useAsDefault = function(x, ..., na.rm = FALSE) { base::range(x, ..., na.rm = na.rm) }, group = "Summary") setGeneric("prod", function(x, ..., na.rm = FALSE) { standardGeneric("prod") }, useAsDefault = function(x, ..., na.rm = FALSE) { base::prod(x, ..., na.rm = na.rm) }, group = "Summary") setGeneric("sum", function(x, ..., na.rm = FALSE) { standardGeneric("sum") }, useAsDefault = function(x, ..., na.rm = FALSE) { base::sum(x, ..., na.rm = na.rm) }, group = "Summary") } @ Now we need some more user-unfriendly helper functions: <>= .Brob.max <- function(x, ..., na.rm=FALSE){ p <- x@positive val <- x@x if(any(p)){ return(brob(max(val[p]))) } else { return(brob(min(val),FALSE)) } } .Brob.prod <- function(x){ p <- x@positive val <- x@x return(brob(sum(val),(sum(p)%%2)==0)) } .Brob.sum <- function(x){ .Brob.sum.allpositive( x[x>0]) - .Brob.sum.allpositive(-x[x<0]) } .Brob.sum.allpositive <- function(x){ if(length(x)<1){return(as.brob(0))} val <- x@x p <- x@positive mv <- max(val) return(brob(mv + log1p(sum(exp(val[-which.max(val)]-mv))),TRUE)) } @ Note the final function that sums its arguments, which are all assumed to be positive, in an intelligent, accurate, and efficient manner. No checking is done (this is not a user-friendly function!). We can define {\tt .Brob.sum()} in terms of this: return the difference between the sum of the positive arguments and and the sum of minus the negative arguments. <>= <>= setMethod("Summary", "brob", function(x, ..., na.rm=FALSE){ switch(.Generic, max = .Brob.max( x, ..., na.rm=na.rm), min = -.Brob.max(-x, ..., na.rm=na.rm), range = cbrob(min(x,na.rm=na.rm),max(x,na.rm=na.rm)), prod = .Brob.prod(x), sum = .Brob.sum(x), stop(paste(.Generic, "not allowed on Brobdingnagian numbers")) ) } ) @ Better check: <>= sum(as.brob(1:100)) - 5050 @ showing acceptable accuracy. \section{Examples of the package in use} We can try to evaluate a factorial. Stirling's approximation is~$n!\sim\sqrt{2\pi n}\,e^{-n}n^n$: <>= stirling <- function(x){sqrt(2*pi*x)*exp(-x)*x^x} @ And this should work seamlessly with Brobs: <>= stirling(100) stirling(as.brob(100)) @ And are they the same? <>= as.numeric(stirling(100)/stirling(as.brob(100))) @ \ldots pretty near. But the great advantage of Brobdingnagian numbers is that they can handle numbers larger than the IEEE limit: <>= stirling(1000) stirling(as.brob(1000)) @ \noindent and this is accurate to about~12 sig figs, which is accurate enough for many purposes. The number of sig figs decreases with progressively larger numbers, essentially because an increasing amount of floating point accuracy is gobbled up by storing the exponent of a large number, and less is left for the mantissa. \section*{Acknowledgments} I gratefully acknowledge the help given to me by many members of the \proglang{R}-help and \proglang{R}-devel lists, especially Martin Maechler, Brian Ripley, and John Chambers. \nocite{*} \bibliography{brob} \end{document}