R/array.R

In Rmpfr: R MPFR - Multiple Precision Floating-Point Reliable

## From an "mpfr" object make an mpfr(Array|Matrix) :

setMethod("dim", "mpfrArray", function(x) x@Dim)
setMethod("dimnames", "mpfrArray", function(x) x@Dimnames)

## 2 basic methods to construct "mpfr - arrays" ( mpfrArray | mpfrMatrix ) :

##' "mpfr" --> "mpfrArray"  --- basically  dim(<mpfr>) <- dd
mpfr2array <- function(x, dim, dimnames=NULL, check=FALSE) {
    if(check) stopifnot(extends((clx <- class(x)), "mpfr"))
    if(is.numeric(dim) && all(dim == (iv <- as.integer(dim)))) {
	rnk <- length(iv)
	if(check) {
	    cl <- if(rnk == 2) "mpfrMatrix" else "mpfrArray"
	    if(extends(clx, "mpfrArray")) x <- as(x, "mpfr")# drop 'Dim', 'Dimnames'
	    if(is.null(dimnames))
		new(cl, x, Dim = iv)
	    else new(cl, x, Dim = iv, Dimnames = dimnames)
	} else { ## faster, non-checking
	    r <- setDataPart(new(if(rnk == 2) "mpfrMatrix" else "mpfrArray"),
			     x, check=FALSE)
	    r@Dim <- iv
	    if(!is.null(dimnames))
		r@Dimnames <- dimnames
##TODO R >= 2.13.2:
##TODO	    else if(.hasSlot(x, "Dimnames")) # has "wrong' Dimnames
	    else if(is(x, "mpfrArray")) # has "wrong' Dimnames
		r@Dimnames <- rep.int(list(NULL), rnk)
	    r
	}
    }
    else if(is.null(dim))
	as.vector(x)
    else
	stop("invalid dimension specified")
}
setMethod("dim<-", signature(x = "mpfr", value = "ANY"),
	  function(x, value) mpfr2array(x, value))


mpfrArray <- function(x, precBits, dim = length(x), dimnames = NULL,
		      rnd.mode = c('N','D','U','Z','A'))
{
    if(!is.atomic(x))
	stop("'x' must be (coercable to) a numeric vector, possibly consider mpfr2array()")
    dim <- as.integer(dim)
    rnd.mode <- toupper(rnd.mode)
    rnd.mode <- match.arg(rnd.mode)

    ml <- .Call(d2mpfr1_list, x, precBits, rnd.mode)
    vl <- prod(dim)
    if (length(x) != vl) {
	if (vl > .Machine$integer.max)
	    stop("'dim' specifies too large an array")
	ml <- rep(ml, length.out = vl)
    }
    new(if(length(dim) == 2) "mpfrMatrix" else "mpfrArray",
	ml, Dim = dim,
	Dimnames = if(is.null(dimnames)) vector("list", length(dim))
		   else dimnames)
}

setAs("array", "mpfr",        function(from) mpfr(from, 128L))
setAs("array", "mpfrArray",   function(from) mpfr(from, 128L))
setAs("matrix", "mpfrMatrix", function(from) mpfr(from, 128L))

## and for "base" functions to work:
as.array.mpfr <- function(x, ...) {
    if(is(x, "mpfrArray")) x else ## is(x, "mpfr") :
    as.array.default(x, ...)
}
as.matrix.mpfr <- function(x, ...) {
    if(is(x, "mpfrMatrix")) x else ## is(x, "mpfr") :
    as.matrix.default(x, ...)
}

## matrix is S3 generic from 'gmp' anyway:
matrix.mpfr <- function (data = NA, nrow = 1, ncol = 1, byrow = FALSE, ...) {
    dim(data) <- c(nrow, ncol)
    if(length(dots <- list(...))) {
	if(!is.null(dn <- dots$dimnames)) {
	    dimnames(data) <- dn # assign and delete from "dots":
	    dots$dimnames <- NULL
	}
	if(nx <- length(dots)) # a simplified  Matrix:::chk.s()
	    warning(sprintf(ngettext(nx,
				     "extra argument %s will be disregarded",
				     "extra arguments %s will be disregarded"),
			sub(")$", '', sub("^list\\(", '', deparse(dots, control=c())))))
    }
    data
}


setMethod("dimnames<-", signature(x = "mpfrArray", value = "ANY"),
	  function(x, value) {
	      if(!is.list(value)) stop("non-list RHS")
	      if(length(value) != length(x@Dim))
		  stop("RHS (new dimnames) differs in length from dim(.)")
	      x@Dimnames <- value
	      x
	  })

setMethod("t", "mpfrMatrix", t.mpfrMatrix <-
	  function(x) {
	      d <- x@Dim; n <- d[1]; m <- d[2]
	      ## These are the indices to get the transpose of m {n x m} :
	      ## ind.t <- function(n,m)rep.int(1:n, rep(m,n)) + n*(0:(m-1))
	      x@Dim <- c(m,n)
	      x@Dimnames <- x@Dimnames[2:1]
	      ## faster than { x@.Data <- x@.Data[rep.int(1:n, rep(m,n)) + n*(0:(m-1))] ; x } :
	      setDataPart(x, getD(x)[rep.int(1:n, rep(m,n)) + n*(0:(m-1))], check=FALSE)
	  })
setMethod("t", "mpfr", t.mpfr <-
	  function(x) { # t(<n-vector>) |-->  {1 x n} matrix
	      r <- new("mpfrMatrix")
	      r@Dim <- c(1L, length(x))
	      ## faster than  { r@.Data <- x@.Data ; r } :
	      setDataPart(r, getD(x), check=FALSE)
	  })

setMethod("aperm", signature(a="mpfrArray"), aperm.mpfrArray <-
	  function(a, perm, resize=TRUE, ...) {
	      stopifnot(1 <= (k <- length(d <- a@Dim)))
	      if(missing(perm)) perm <- k:1
	      else stopifnot(length(perm <- as.integer(perm)) == k, 1 <= perm, perm <= k)
	      if(!resize)
		  stop("'resize != TRUE is not (yet) implemented for 'mpfrArray'")
	      a@Dim <- d[perm]
	      a@Dimnames <- a@Dimnames[perm]
	      ii <- c(aperm(array(1:prod(d), dim=d), perm=perm, resize=FALSE))
	      ## faster than  { a@.Data <- a@.Data[ ii ] ; a } :
	      setDataPart(a, getD(a)[ ii ], check=FALSE)
	  })


## `` drop the  dim() part '' :
setMethod("as.vector", "mpfrArray", function(x) as(x, "mpfr"))
## a "vector" in  *one* sense at least, and "mpfr" does extend "vector":
setAs("mpfrArray", "vector", function(from) as(from, "mpfr"))

.toNum <- function(from, rnd.mode) { ## <- must have only 'from'
    if(is.null(dn <- dimnames(from)) || identical(dn, list(NULL,NULL)))
	## --> result has NULL dimnames
	structure(.Call(mpfr2d, from, rnd.mode), dim = dim(from))
    else
	structure(.Call(mpfr2d, from, rnd.mode), dim = dim(from), dimnames = dn)
}
## to be used in setAs(),  must have only 'from' argument:
.toNum1 <- function(from) .toNum(from, rnd.mode="N")

toNum <- function(from, rnd.mode = c('N','D','U','Z','A')) {
    stopifnot(is.character(rnd.mode <- toupper(rnd.mode)))
    rnd.mode <- match.arg(rnd.mode)
    .toNum(from, rnd.mode)
}

setAs("mpfrArray", "array", .toNum1)

setAs("mpfrMatrix", "matrix", .toNum1)

setAs("mpfrArray", "matrix", function(from) {
    if(length(dim(from)) != 2)
	stop("dim(.) != 2  ==> cannot be coerced to 'matrix'")
    toNum(from)
})

print.mpfrArray <-
    function(x, digits = NULL, drop0trailing = FALSE, right = TRUE,
             ##				      -----
             ## would like    'drop0... = TRUE', but that's only ok once we have a
             ## format() allowing to "jointly format a column"
             max.digits = getOption("Rmpfr.print.max.digits", 999L),
             exponent.plus = getOption("Rmpfr.print.exponent.plus", TRUE),
             ...)
{
    stopifnot(is(x, "mpfrArray"), is.null(digits) || digits >= 2)
    ## digits = NULL --> the inherent precision of x will be used
    n <- length(x)
    ch.prec <-
	if(n >= 1) {
	    rpr <- range(.getPrec(x))
	    paste("of precision ", rpr[1],
		   if(rpr[1] != rpr[2]) paste("..",rpr[2]), " bits")
	}
    cl <- class(x)
    p0 <- function(...) paste(..., sep="")
    cat(p0("'",cl,"'"), "of dim(.) = ",
	p0("(",paste(x@Dim, collapse=", "),")"),
	ch.prec, "\n")
    if(n >= 1) {
	## FIXME: really need a 'format' method for mpfrArrays
	## -----  which properly aligns columns !!

	## Build character array fx, and print that

        ## drop arguments for print.default(*):
	lFormat <- function(x, na.print, print.gap, max, useSource, ...)
	    formatMpfr(x, digits=digits, max.digits=max.digits,
                   drop0trailing=drop0trailing, exponent.plus=exponent.plus,
		   ...)
	fx <- lFormat(x, ...)
	dim(fx) <- dim(x)
	dimnames(fx) <- dimnames(x)
	print(fx, ..., right=right, quote = FALSE)
    }
    invisible(x)
}
setMethod(show, "mpfrArray", function(object) print.mpfrArray(object))

## FIXME : should happen in C, where we could "cut & paste" much of
## -----  do_matprod() and matprod() from ~/R/D/r-devel/R/src/main/array.c
##/* "%*%" (op = 0), crossprod (op = 1) or tcrossprod (op = 2) */
.matmult.R <- function(x,y, op = 0L, fPrec = 1, precBits = fPrec * max(getPrec(x), getPrec(y)))
{
    if(!(is.numeric(x) || is.mpfr(x)))
	stop("'x' must be numeric or mpfr(Matrix)")
    sym <- missing(y)
    if (sym && (op > 0L)) y <- x
    else if(!(is.numeric(y) || is.mpfr(y)))
	stop("'y' must be numeric or mpfr(Matrix)")
    ldx <- length(dx <- dim(x))
    ldy <- length(dy <- dim(y))
    ## "copy, paste & modify" from  do_matprod():
    if (ldx != 2 && ldy != 2) {		#* x and y non-matrices */
	if (op == 0L) {
	    nrx <- 1L; ncx <- length(x)
	} else {
	    nrx <- length(x); ncx <- 1L
	}
	nry <- length(y)
	ncy <- 1L
    }
    else if (ldx != 2) {		#* x not a matrix */
	nry <- dy[1]
	ncy <- dy[2]
	nrx <- ncx <- 0L
	if (op == 0L) {
	    if (length(x) == nry) {	#* x as row vector */
		nrx <- 1L
		ncx <- nry # == length(x)
	    }
	    else if (nry == 1) {	#* x as col vector */
		nrx <- length(x)
		ncx <- 1L # == nry
	    }
	}
	else if (op == 1L) { #* crossprod
	    if (length(x) == nry) {	#* x is a col vector */
		nrx <- nry # = length(x)
		ncx <- 1L
	    }
	}
	else { # op == 2L: tcrossprod
	    if (length(x) == ncy) {	#* x as row vector */
		nrx <- 1L
		ncx <- ncy # == length(x)
	    }
	    else if (ncy == 1) {	#* x as col vector */
		nrx <- length(x)
		ncx <- 1L # == ncy
	    }
	}
    }
    else if (ldy != 2) {		#* y not a matrix */
	nrx <- dx[1]
	ncx <- dx[2]
	nry <- ncy <- 0L
	if (op == 0L) {
	    if (length(y) == ncx) {	#* y as col vector */
		nry <- ncx # = length(y)
		ncy <- 1L
	    }
	    else if (ncx == 1) {	#* y as row vector */
		nry <- 1L # = ncx
		ncy <- length(y)
	    }
	}
	else if (op == 1L) { #* crossprod
	    if (length(y) == nrx) {	#* y is a col vector */
		nry <- nrx # = length(y)
		ncy <- 1L
	    }
	}
	else { # op == 2L: tcrossprod	y is a col vector
	    nry <- length(y)
	    ncy <- 1L
	}
    }
    else {				#* x and y matrices */
	nrx <- dx[1]
	ncx <- dx[2]
	nry <- dy[1]
	ncy <- dy[2]
    }
    ##* nr[ow](.) and nc[ol](.) are now defined for x and y */

    z <- new("mpfrMatrix")
    z0 <- as(0, "mpfr")

    if (op == 0L) { ## %*%
	if (ncx != nry) stop("non-conformable arguments")

	z@Dim <- c(nrx, ncy)
	z@.Data <- vector("list", nrx*ncy)
	if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
	    j <- 0L:(ncx - 1L)
	    for(i in 1:nrx) {
		for (k in 0L:(ncy - 1L))
		    z[i + k * nrx] <-
			## sum(x[i + j * nrx] * y[1L+ j + k * nry])
			new("mpfr",
			    .Call(R_mpfr_sumprod, # --> ../src/Summary.c
				  x[i + j * nrx], y[1L+ j + k * nry],
				  precBits, alternating=FALSE))
	    }
	}
	else	    #/* zero-extent operations should return zeroes */
	    for(i in seq_len(nrx*ncy)) z[i] <- z0
    }
    else if (op == 1L) { ## crossprod() :  x' %*% y
	if (nrx != nry) stop("non-conformable arguments")

	z@Dim <- c(ncx, ncy)
	z@.Data <- vector("list", ncx*ncy)
	if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
	    j <- 1L:nrx
	    for(i in 0L:(ncx - 1L)) {
		for (k in 0L:(ncy - 1L))
		    z[1L +i + k * ncx] <-
			## sum(x[j + i * nrx] * y[j + k * nry])
			new("mpfr",
			    .Call(R_mpfr_sumprod,
				  x[j + i * nrx], y[j + k * nry],
				  precBits, alternating=FALSE))

	    }
	} else
	    for(i in seq_len(ncx*ncy)) z[i] <- z0
    }
    else { ## op == 2L :	 tcrossprod() :	 x %*% y'
	if (ncx != ncy) stop("non-conformable arguments")

	z@Dim <- c(nrx, nry)
	z@.Data <- vector("list", nrx*nry)
	if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0)
	    for(i in seq_len(nrx)) {
		j <- 0L:(ncx - 1L)
		for (k in 0L:(nry - 1L))
		    z[i + k * nrx] <-
			## sum(x[i + j * nrx] * y[1L +k + j * nry])
			new("mpfr",
			    .Call(R_mpfr_sumprod,
				  x[i + j * nrx], y[1L +k + j * nry],
				  precBits, alternating=FALSE))

	    }
	else
	    for(i in seq_len(nrx*nry)) z[i] <- z0
    }
    z
} ## .matmult.R()

## "FIXME"?  make working also with "Matrix" class matrices ..
##                             ----------------------------

## An 'explicit' %*% function (with >= 2 arguments) :
matmult <- function(x,y, ...) .matmult.R(x,y, op = 0L, ...)

setMethod("%*%", signature(x = "mpfrMatrix", y = "mpfrMatrix"),
	  function(x,y) .matmult.R(x,y, op= 0L))
setMethod("%*%", signature(x = "mpfrMatrix", y = "mpfr"),
	  function(x,y) .matmult.R(x,y, op= 0L))
setMethod("%*%", signature(x = "mpfr", y = "mpfrMatrix"),
	  function(x,y) .matmult.R(x,y, op= 0L))
setMethod("%*%", signature(x = "mpfr", y = "mpfr"),
	  function(x,y) .matmult.R(x,y, op= 0L))
## These cover vectors, etc (!) :
setMethod("%*%", signature(x = "mpfr", y = "Mnumber"),
	  function(x,y) .matmult.R(x,y, op= 0L))
setMethod("%*%", signature(x = "Mnumber", y = "mpfr"),
	  function(x,y) .matmult.R(x,y, op= 0L))


setMethod("crossprod", signature(x = "mpfrMatrix", y = "mpfrMatrix"),
	  function(x,y, ...) .matmult.R(x,y, op= 1L, ...))
setMethod("crossprod", signature(x = "mpfrMatrix", y = "mpfr"),
	  function(x,y, ...) .matmult.R(x,y, op= 1L, ...))
setMethod("crossprod", signature(x = "mpfr", y = "mpfrMatrix"),
	  function(x,y, ...) .matmult.R(x,y, op= 1L, ...))
setMethod("crossprod", signature(x = "mpfr", y = "mpfr"),
	  function(x,y, ...) .matmult.R(x,y, op= 1L, ...))
setMethod("crossprod", signature(x = "mpfr", y = "Mnumber"),
	  function(x,y, ...) .matmult.R(x,y, op= 1L, ...))
setMethod("crossprod", signature(x = "Mnumber", y = "mpfr"),
	  function(x,y, ...) .matmult.R(x,y, op= 1L, ...))
## one argument-case: [FIXME: not copying 'x' and using x_i^2  is more efficient]
setMethod("crossprod", signature(x = "mpfr", y = "missing"),
	  function(x,y, ...) .matmult.R(x,x, op= 1L, ...))

setMethod("tcrossprod", signature(x = "mpfrMatrix", y = "mpfrMatrix"),
	  function(x,y, ...) .matmult.R(x,y, op= 2L, ...))
setMethod("tcrossprod", signature(x = "mpfrMatrix", y = "mpfr"),
	  function(x,y, ...) .matmult.R(x,y, op= 2L, ...))
setMethod("tcrossprod", signature(x = "mpfr", y = "mpfrMatrix"),
	  function(x,y, ...) .matmult.R(x,y, op= 2L, ...))
setMethod("tcrossprod", signature(x = "mpfr", y = "mpfr"),
	  function(x,y, ...) .matmult.R(x,y, op= 2L, ...))
setMethod("tcrossprod", signature(x = "mpfr", y = "Mnumber"),
	  function(x,y, ...) .matmult.R(x,y, op= 2L, ...))
setMethod("tcrossprod", signature(x = "Mnumber", y = "mpfr"),
	  function(x,y, ...) .matmult.R(x,y, op= 2L, ...))
## one argument-case: [FIXME: not copying 'x' and using x_i^2  is more efficient]
setMethod("tcrossprod", signature(x = "mpfr", y = "missing"),
	  function(x,y, ...) .matmult.R(x,x, op= 2L, ...))


.mpfrA.subset <- function(x,i,j, ..., drop) {
    nA <- nargs()
    if(getOption("verbose"))
	message(sprintf("nargs() == %d  mpfrArray indexing ... ", nA))

    r <- getD(x) # the data part, a list()
    if(nA == 2) ## A[i]
	return(new("mpfr", r[i]))
    ## else: nA != 2 : nA > 2 -
    dim(r) <- dim(x)
    dimnames(r) <- dimnames(x)
    r <- r[i,j, ..., drop=drop]
    if(drop && is.null(dim(r)))
	new("mpfr", r)
    else {
	D <- if(is.null(dr <- dim(r))) # ==> drop is FALSE; can this happen?
	    rep.int(1L, length(r)) else dr
	x@Dim <- D
	x@Dimnames <- if(is.null(dn <- dimnames(r)))
	    vector("list", length(D)) else dn
	if(length(D) == 2 && !inherits(x, "mpfrMatrix"))
	    ## low-level "coercion" from mpfrArray to *Matrix :
	    attr(x,"class") <- getClass("mpfrMatrix")@className
	attributes(r) <- NULL
	setDataPart(x, r, check=FALSE)
    }
}

## "["
setMethod("[", signature(x = "mpfrArray", i = "ANY", j = "ANY", drop = "ANY"),
	  .mpfrA.subset)

## this signature needs a method here, or it triggers the one for "mpfr"
setMethod("[", signature(x = "mpfrArray", i = "ANY", j = "missing",
			 drop = "missing"),
	  .mpfrA.subset)


.mA.subAssign <- function(x,i,j,..., value, n.a, isMpfr)
{
    ## n.a :=== nargs() -- in the calling "[<-" method --
    r <- getD(x)
    if(n.a >= 4) {
	## A[i,j]  /  A[i,]  /	A[,j]	but not A[i]
	## A[i,j,k] <- v : n.a == 5
	dim(r) <- dim(x)
	dimnames(r) <- dimnames(x)
	if(!isMpfr)
	   value <- mpfr(value, precBits =
			 pmax(getPrec(value),
			      .getPrec(if(n.a == 4) r[i,j] else r[i,j, ...]))
			 )
	vD <- getD(value)
	if(n.a == 4) {
	    r[i,j] <- vD
	} else { ## n.a >= 5
	    r[i, j, ...] <- vD
	}
	attributes(r) <- NULL
    }
    else if(n.a %in% c(2,3)) { ##  A [ i ] <- v  //   A[] <- v
	if(!isMpfr)
	    value <- mpfr(value, precBits = pmax(getPrec(value), .getPrec(r[i])))
        if(n.a == 3L)
            r[i] <- value
        else ## n.a == 2:
            r[] <- value
    }
    else { ## n.a <= 1
	stop(sprintf("nargs() == %d  mpfrArray[i,j] <- value  __ SHOULD NOT HAPPEN!",
                     n.a))
    }
    setDataPart(x, r, check=FALSE)
}## .mA.subAssign

## "[<-" :
## -------
## E.g., for A[1,,2] <- V
## these are to trigger before the  ("mpfr", i,j, "mpfr")  [ ./mpfr.R ] does
for(it in c("ANY", "missing"))
    for(jt in c("ANY", "missing"))
    setReplaceMethod("[", signature(x = "mpfrArray", i = it, j = jt, value = "mpfr"),
		     function(x,i,j,..., value)
		     .mA.subAssign(x,i=i,j=j,...,value=value,
				   n.a=nargs(), isMpfr = TRUE))
## non-"mpfr" value
for(it in c("ANY", "missing"))
    for(jt in c("ANY", "missing"))
    setReplaceMethod("[", signature(x = "mpfrArray", i = it, j = jt, value = "ANY"),
		     function(x,i,j, ..., value)
		     .mA.subAssign(x,i=i,j=j,...,value=value,
				   n.a=nargs(), isMpfr = FALSE))
rm(it,jt)

## In the Matrix package we have  Diagonal()  for *constructing* a diagonalMatrix;
## in any case, we do only want to support the   diag(<matrix>)    case.

setMethod("diag", signature(x = "mpfrMatrix"),
	  function(x, nrow, ncol) {
	      n <- min(dim(x)); i <- seq_len(n); x[cbind(i,i)] })

setMethod("diag<-", signature(x = "mpfrMatrix"),
	  function(x, value) {
	      n <- min(dim(x)); i <- seq_len(n); x[cbind(i,i)] <-  value; x })



###-----------

setGeneric("cbind", signature = "...")# -> message about override & deparse.level
setGeneric("rbind", signature = "...")

## inside such cbind() / rbind() S4 methods,  match.call() does *not* work correctly,
## this works *only* for top-level calls :
bind_match.call <- function() sys.call(1L)
## so use our "hack" :
bind_match.call <- function() {
    nc <- length(scs <- sys.calls()) # last one is bind_match.call() itself
    ## want the one call *above*  standardGeneric("...") :
    if(is.symbol(fn <- scs[[nc-1L]][[1L]])) { # e.g. 'cbind'
        Gcall <- call("standardGeneric", as.character(fn)) # e.g. standardGeneric("cbind")
        i. <- which(vapply(scs, identical, NA, Gcall))
        scs[[if(!length(i.) || i. < 2L) 1L else i. - 1L ]]
    } else # try "better"
        match.call()
}


setMethod("cbind", "Mnumber",
	  function(..., deparse.level = 1) {
	      args <- list(...)
	      if(all(vapply(args, is.atomic, NA)))
		  return( base::cbind(..., deparse.level = deparse.level) )
	      ## else: at least one is "mpfr(Matrix/Array)"

	      if(any(vapply(args, is.character, NA))) {
		  ## result will be  <character> matrix !
		  isM <- vapply(args, is, NA, class2 = "mpfr")
		  args[isM] <- lapply(args[isM], as, Class = "character")
		  return(do.call(base::cbind,
				 c(args, list(deparse.level=deparse.level))))

	      } else if(any(vapply(args, is.complex, NA))) {
		  ## result will be  <complex> matrix;
		  ## in the future <complex_mpfr>  ???

		  stop("cbind(...) of 'complex' and 'mpfr' objects is not implemented")
		  ## give at least warning !!
	      }
	      ## else
	      L <- function(a) if(is.numeric(n <- nrow(a))) n else length(a)
	      W <- function(a) if(is.numeric(n <- ncol(a))) n else 1L
	      ## the number of rows of the result : {for now require integer}
	      NR <- max(lengths <- vapply(args, L, integer(1)))
	      NC <- sum(widths	<- vapply(args, W, integer(1)))
	      r <- setDataPart(new("mpfrMatrix"), vector("list", NR*NC))
	      r@Dim <- as.integer(c(NR, NC))
	      hasDim <- !vapply(args, function(a) is.null(dim(a)), NA)
	      do.colnames <- deparse.level || any(hasDim)
	      if(do.colnames) {
		  nms <- character(NC)
		      ## help(cbind)  has in 'Value' :
		      ##      For ‘cbind’ (‘rbind’) the column (row) names are taken from the
		      ##      colnames (rownames) of the arguments if these are matrix-like.
		      ##      Otherwise from the names of the arguments or where those are not
		      ##      supplied and ‘deparse.level > 0’, by deparsing the expressions
		      ##      given, for ‘deparse.level = 1’ only if that gives a sensible name
		      ##      (a symbol, see is.symbol).
		  nV <- names(widths) # == names(args), possibly NULL
		  hasV <- !is.null(nV)
		  ## argN <- substitute(...)## "fails" here same as match.call()
		  fcall <- bind_match.call()
		  if(!missing(deparse.level)) # must remove: it could be "anywhere"
		      fcall <- fcall[names(fcall) != "deparse.level"]
		  ## cat("fcall: "); str(fcall)
		  ## browser()
		  ## argN <- fcall[-1] # but that makes 1st arg into fn.name!
		  ## vapply(fcall[-1], deparse1, "") ## is what we'd need, incl. *names*
		  ## not ok when called as  selectMethod("cbind","mpfr")(x, ....)
                  fcall.ok <- (length(fcall) == 1L + ...length()) ## == 1 + length(args)
	      }
	      j <- 0
	      prec <- .Machine$double.digits
	      for(ia in seq_along(args)) {
		  w <- widths[ia]
		  a <- args[[ia]]
                  isM <- hasDim[[ia]] # == !is.null(dim(a)) ; true iff matrix-like
		  if(is.mpfr(a)) {
		      prec <- max(prec, .getPrec(a))
		  } else { ## not "mpfr"
		      a <- mpfr(a, prec)
		  }
		  if((li <- lengths[ia]) != 1 && li != NR) { ## recycle
		      if(isM)
			  stop("number of rows of matrices must match")
		      ## else
		      if(NR %% li)
			  warning("number of rows of result is not a multiple of vector length")
		      a <- a[rep(seq_len(li), length.out = NR)]
		  }
		  ii <- j+ seq_len(w)
		  r[, ii] <- a
		  if(do.colnames) {
		      nms[ii] <-
			  if(isM)
			      colnames(a) %||% ""
			  else {
			      if(hasV && nzchar(n. <- nV[[ia]]))
				  n.
			      else if(fcall.ok) { ## deparsed argument "arg"[[ia]]
				  a <- fcall[[ia+1L]]
				  if(is.symbol(a) || deparse.level == 2)
				      deparse1(a)
				  else ""
			      } else ""
			  }
		  }
		  j <- j + w
	      }
	      if(do.colnames && any(nzchar(nms)))
		  r@Dimnames[[2L]] <- nms
	      r
	  })

setMethod("rbind", "Mnumber",
	  function(..., deparse.level = 1) {
	      args <- list(...)
	      if(all(vapply(args, is.atomic, NA)))
		  return( base::rbind(..., deparse.level = deparse.level) )
	      ## else: at least one is "mpfr(Matrix/Array)"

	      if(any(vapply(args, is.character, NA))) {
		  ## result will be  <character> matrix !
		  isM <- vapply(args, is, NA, class2 = "mpfr")
		  args[isM] <- lapply(args[isM], as, Class = "character")
		  return(do.call(base::rbind,
				 c(args, list(deparse.level=deparse.level))))

	      } else if(any(vapply(args, is.complex, NA))) {
		  ## result will be  <complex> matrix;
		  ## in the future <complex_mpfr>  ???

		  stop("rbind(...) of 'complex' and 'mpfr' objects is not implemented")
		  ## give at least warning !!
	      }
	      ## else
	      L <- function(a) if(is.numeric(n <- nrow(a))) n else 1L
	      W <- function(a) if(is.numeric(n <- ncol(a))) n else length(a)
	      ## the number of rows of the result : {for now require integer}
	      NR <- sum(lengths <- vapply(args, L, integer(1)))
	      NC <- max(widths	<- vapply(args, W, integer(1)))
	      r <- setDataPart(new("mpfrMatrix"), vector("list", NR*NC))
	      r@Dim <- as.integer(c(NR, NC))
	      hasDim <- !vapply(args, function(a) is.null(dim(a)), NA)
	      do.rownames <- deparse.level || any(hasDim)
	      if(do.rownames) {
		  nms <- character(NR)
		  nV <- names(lengths) # == names(args), possibly NULL
		  hasV <- !is.null(nV)
		  fcall <- bind_match.call()
		  if(!missing(deparse.level)) # must remove: it could be "anywhere"
		      fcall <- fcall[names(fcall) != "deparse.level"]
		  ## not ok when called as  selectMethod("cbind","mpfr")(x, ....)
                  fcall.ok <- (length(fcall) == 1L + ...length()) ## == 1 + length(args)
	      }
	      i <- 0
	      prec <- .Machine$double.digits
	      for(ia in seq_along(args)) {
		  le <- lengths[ia]
		  a <- args[[ia]]
		  isM <- hasDim[[ia]] # == !is.null(dim(a)) ; true iff matrix-like
		  if(is.mpfr(a)) {
		      prec <- max(prec, .getPrec(a))
		  } else { ## not "mpfr"
		      a <- mpfr(a, prec)
		  }
		  if((wi <- widths[ia]) != 1 && wi != NC) { ## recycle
		      if(isM)
			  stop("number of rows of matrices must match")
		      ## else
		      if(NC %% wi)
			  warning("number of columns of result is not a multiple of vector length")
		      a <- a[rep(seq_len(wi), length.out = NC)]
		  }
		  ii <- i+ seq_len(le)
		  r[ii, ] <- a
		  if(do.rownames) {
		      nms[ii] <-
			  if(isM)
			      rownames(a) %||% ""
			  else {
			      if(hasV && nzchar(n. <- nV[[ia]]))
				  n.
			      else if(fcall.ok) { ## deparsed argument "arg"[[ia]]
				  a <- fcall[[ia+1L]]
				  if(is.symbol(a) || deparse.level == 2)
				      deparse1(a)
				  else ""
			      } else ""
			  }
		  }
		  i <- i + le
	      }
	      if(do.rownames && any(nzchar(nms)))
		  r@Dimnames[[1L]] <- nms
	      r
	  })

unlistMpfr <- function(x, recursive = FALSE, use.names = TRUE)  {
    ## an "unlist(.)" for mpfr contents:
    if(recursive) stop("'recursive = TRUE' is not implemented (yet).")
    n <- sum(lengths(x))
    ans <- mpfr(numeric(n), precBits=3L)# dummy to fill
    ans@.Data <- unlist(lapply(x, slot, ".Data"), use.names=use.names)
    ans
}


##-- Original in ~/R/D/r-devel/R/src/library/base/R/apply.R :
##
applyMpfr <- function(X, MARGIN, FUN, ...)
{
    FUN <- match.fun(FUN)

    ## Ensure that X is an array object
    dl <- length(dim(X))
    if(!dl) stop("dim(X) must have a positive length")
##-     if(is.object(X))
##-	X <- if(dl == 2L) as.matrix(X) else as.array(X)
    ## now record dim as coercion can change it
    ## (e.g. when a data frame contains a matrix).
    d <- dim(X)
    dn <- dimnames(X)
    ds <- seq_len(dl)

    ## Extract the margins and associated dimnames

    if (is.character(MARGIN)) {
	if(is.null(dnn <- names(dn))) # names(NULL) is NULL
	   stop("'X' must have named dimnames")
	MARGIN <- match(MARGIN, dnn)
	if (any(is.na(MARGIN)))
	    stop("not all elements of 'MARGIN' are names of dimensions")
    }
    s.call <- ds[-MARGIN]
    s.ans  <- ds[MARGIN]
    d.call <- d[-MARGIN]
    d.ans  <- d[MARGIN]
    dn.call<- dn[-MARGIN]
    dn.ans <- dn[MARGIN]
    ## dimnames(X) <- NULL

    array <- function(data, dim = length(data), dimnames = NULL) {
	dim(data) <- dim
	if(!is.null(dimnames)) dimnames(data) <- dimnames
	data
    }

    ## do the calls

    d2 <- prod(d.ans)
    if(d2 == 0L) {
	## arrays with some 0 extents: return ``empty result'' trying
	## to use proper mode and dimension:
	## The following is still a bit `hackish': use non-empty X
	newX <- array(vector(typeof(X), 1L), dim = c(prod(d.call), 1L))
	ans <- forceAndCall(1, FUN, if(length(d.call) < 2L) newX[,1] else
		   array(newX[, 1L], d.call, dn.call), ...)
	return(if(is.null(ans)) ans else if(length(d.ans) < 2L) ans[1L][-1L]
	       else array(ans, d.ans, dn.ans))
    }
    ## else
    newX <- aperm(X, c(s.call, s.ans))
    dim(newX) <- c(prod(d.call), d2)
    ans <- vector("list", d2)
    if(length(d.call) < 2L) {# vector
	if (length(dn.call)) dimnames(newX) <- c(dn.call, list(NULL))
	for(i in 1L:d2) {
	    tmp <- forceAndCall(1, FUN, newX[,i], ...)
	    if(!is.null(tmp)) ans[[i]] <- tmp
	}
    } else
       for(i in 1L:d2) {
	   tmp <- forceAndCall(1, FUN, array(newX[,i], d.call, dn.call), ...)
	   if(!is.null(tmp)) ans[[i]] <- tmp
	}

    ## answer dims and dimnames

    ans.list <- !is(ans[[1L]], "mpfr") ##- is.recursive(ans[[1L]])
    l.ans <- length(ans[[1L]])

    ans.names <- names(ans[[1L]])
    if(!ans.list)
	ans.list <- any(lengths(ans) != l.ans)
    if(!ans.list && length(ans.names)) {
	all.same <- vapply(ans, function(x) identical(names(x), ans.names), NA)
	if (!all(all.same)) ans.names <- NULL
    }
    len.a <- if(ans.list) d2 else length(ans <- unlistMpfr(ans))
    if(length(MARGIN) == 1L && len.a == d2) {
	names(ans) <- if(length(dn.ans[[1L]])) dn.ans[[1L]] # else NULL
	ans
    }
    else if(len.a == d2)
	array(ans, d.ans, dn.ans)
    else if(len.a && len.a %% d2 == 0L) {
	if(is.null(dn.ans)) dn.ans <- vector(mode="list", length(d.ans))
	dn1 <- list(ans.names)
	if(length(dn.call) && !is.null(n1 <- names(dn <- dn.call[1])) &&
	   nzchar(n1) && length(ans.names) == length(dn[[1]]))
	    names(dn1) <- n1
	dn.ans <- c(dn1, dn.ans)
	array(ans, c(len.a %/% d2, d.ans),
	      if(!is.null(names(dn.ans)) || !all(vapply(dn.ans, is.null, NA)))
		  dn.ans)
    } else
	ans
}

setGeneric("apply")
setMethod ("apply", "mpfrArray", applyMpfr)

setMethod("colSums", "mpfrArray", function(x, na.rm = FALSE, dims = 1, ...) {
    stopifnot((rnk <- length(dim(x))) >= 2, 1 <= dims, dims <= rnk - 1)
    applyMpfr(x, (dims+1):rnk, sum)
})
setMethod("colMeans", "mpfrArray", function(x, na.rm = FALSE, dims = 1, ...) {
    stopifnot((rnk <- length(dim(x))) >= 2, 1 <= dims, dims <= rnk - 1)
    applyMpfr(x, (dims+1):rnk, mean)
})
setMethod("rowSums", "mpfrArray", function(x, na.rm = FALSE, dims = 1, ...) {
    stopifnot((rnk <- length(dim(x))) >= 2, 1 <= dims, dims <= rnk - 1)
    applyMpfr(x, 1:dims, sum)
})
setMethod("rowMeans", "mpfrArray", function(x, na.rm = FALSE, dims = 1, ...) {
    stopifnot((rnk <- length(dim(x))) >= 2, 1 <= dims, dims <= rnk - 1)
    applyMpfr(x, 1:dims, mean)
})

## Cut'n'paste from  ~/R/Pkgs/Matrix/R/Auxiliaries.R {FIXME? load Matrix:::mkDet}
mkDet <- function(d, logarithm = TRUE, ldet = sum(log(abs(d))),
		  sig = -1L+2L*as.integer(prod(sign(d)) >= 0))
{		# sig: -1 or +1 (not 0 !)
    modulus <- if (logarithm) ldet else exp(ldet)
    attr(modulus, "logarithm") <- logarithm
    val <- list(modulus = modulus, sign = sig)
    class(val) <- "det"
    val
}

## S3 method instead of S4, as  base::determinant is S3 generic
determinant.mpfrMatrix <-
    function(x, logarithm = TRUE, asNumeric = (d[1] > 3),
	     precBits = max(.getPrec(x)), ...)
{
    d <- x@Dim
    if(d[1] != d[2]) stop("'x' must ba a square matrix")
    if((n <- d[1]) == 0) determinant(matrix(1,0,0), logarithm=logarithm)
    else if(n == 1)
	mkDet(x[1], logarithm=logarithm)
    else { ## n x n,	for  n >= 2
	if(asNumeric)
	    return(determinant(asNumeric(x), logarithm=logarithm, ...))
	## else use recursive (Care: horribly slow for non-small n!)
	Det <- function(x, n = dim(x)[1]) {
	    if(n == 1) x[1]
	    else if(n == 2) x[1]*x[4] - x[2]*x[3]
	    else {
		a <- mpfr(numeric(n), precBits=3L) # dummy to fill
		n1 <- n-1L
		for(i in seq_len(n)) {
		    a[i] <- Det(x[-i,-1], n=n1)
		}
		## sum(x[,1] * a),  faster :
		new("mpfr",
		    .Call(R_mpfr_sumprod, x[,1], a,
			  precBits, alternating=TRUE))
	    }
	}
	mkDet(Det(x, n=n), logarithm=logarithm)
    }
}

## Only needed for S4 determinant(), not for S3 one:
## The ``Right Thing'' to do :
## base::det() calls [base::]determinant();
## our det() should call our determinant() :
## det <- base::det
## environment(det) <- environment()## == asNamespace("Rmpfr")

if(FALSE) {
## This will become easy, once we have  outer(...)  working, basically almost ==
## base::.kronecker                     ~~~~~~~~~
## ------^---------
setMethod("kronecker", signature(X = "mpfrMatrix", Y = "mpfrMatrix"),
	  function (X, Y, FUN = "*", make.dimnames = FALSE, ...)
      {
	  ydim <- Y@Dim
	  rprec <- max(.getPrec(X),.getPrec(Y))
	  xx <- .......
	  mpfr2array(xx, dim = X@Dim * ydim)
      })
}

scale.mpfrMatrix <- scale.default
## essential, so that colMeans() is using "our" colMeans :
environment(scale.mpfrMatrix) <- environment()# = the "Rmpfr" namespace

### norm() - are "lifted" from ~/R/Pkgs/Matrix/R/sparseMatrix.R :
##  "FIXME": ideally should be part of the setGenericImplicit("norm",..)
setMethod("norm", signature(x = "ANY", type = "missing"),
          function (x, type, ...) norm(x, type = "O", ...))
setMethod("norm", signature(x = "mpfrMatrix", type = "character"),
	  function(x, type, ...) {
	      type <- toupper(substr(type[1], 1, 1))
	      switch(type,  ##  max(<empty>, 0)  |-->  0
		     "O" = ,
                     "1" = max(colSums(abs(x)), 0), ## One-norm (L_1)
		     "I" = max(rowSums(abs(x)), 0), ## L_Infinity
		     "F" = sqrt(sum(x^2)), ## Frobenius
		     "M" = max(abs(x), 0), ## Maximum modulus of all
		     ## otherwise:
		     stop("invalid 'type'"))
	  })


setMethod("head", signature(x = "mpfrMatrix"), utils::head.matrix)
setMethod("tail", signature(x = "mpfrMatrix"), utils::tail.matrix)

## Workaround fact that base::outer() using tcrossprod() does not dispatch (but did on %*% !?!?!?!)
environment(outer) <- environment() # and export and document >>  ../man/base-copies.Rd