Question

I'm mucking with some number theory & would like to obtain a bigz - compatible version of combn (function to calculate comb(n,j) ). bigz is the class name of extended integers the gmp package generates and processes.
I know I can get away with converting bigz to double and working with "integers as floats", at least until I get to some really large integers, but if anyone's written some combinatoric functions compatible with gmp, I'd appreciate a pointer. (otherwise I'll go off and try to write them on my own :-( )

Was it helpful?

Solution

OK, so here's the modded combn code that runs with bigz . Note that input validation is incomplete, and in particular validation that a FUN function has a bigz method is not yet implemented. My final version, which I'll hand off to the maintainer of gmp, (I hope!!) will implement this as a "combin.bigz" method rather than a standalone function.

combnz <- function (x, m, FUN = NULL, simplify = TRUE, ...) 
{
    stopifnot(length(m) == 1L)
    if (m < 0) 
        stop("m < 0", domain = NA)
    # Really not much point in creating a sequence 1:gigundo_bigz
    #if (is.numeric(x) && length(x) == 1L && x > 0 && trunc(x) ==  x) 
        #x <- seq_len(x)
    n <- length(x)
    if (n < m) 
        stop("n < m", domain = NA)
    m <- as.integer(m)
    e <- 0
    h <- m
    a <- seq_len(m)
    nofun <- is.null(FUN)
    if (!nofun && !is.function(FUN)) 
        stop("'FUN' must be a function or NULL")
# TODO: add a " grep('bigz',methods(FUN) )" to verify there's a method.
    len.r <- length(r <- if (nofun) x[a] else FUN(x[a], ...))
    count <- as.integer(round(choose(n, m)))
    if (simplify) {
        dim.use <- if (nofun) 
            c(m, count)
        else {
            d <- dim(r)
            if (length(d) > 1L) 
                c(d, count)
            else if (len.r > 1L) 
                c(len.r, count)
            else c(d, count)
        }
    }
    if (simplify) {
        out <- matrix(r, nrow = len.r, ncol = count)
    }
    else {
        out <- vector("list", count)
        out[[1L]] <- r
    }
    if (m > 0) {
        i <- 2L
        nmmp1 <- n - m + 1L
        while (a[1L] != nmmp1) {
            if (e < n - h) {
                h <- 1L
                e <- a[m]
                j <- 1L
            }
            else {
                e <- a[m - h]
                h <- h + 1L
                j <- 1L:h
            }
            a[m - h + j] <- e + j
            r <- if (nofun) 
                x[a]
            else FUN(x[a], ...)
            if (simplify)
# bigz doesn't  handle replacement the way regular matrices do.  So...
                out[1:len.r,i] <- r
                # out[, i] <- r
# I suspect a list element ain't gonna care
            else out[[i]] <- r
            i <- i + 1L
        }
    }
    if (simplify) 
        array(out, dim.use)
    else out
}
Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top