spectralnorm.R

# ------------------------------------------------------------------
# The Computer Language Shootout
# http://shootout.alioth.debian.org/
#
# Contributed by Leo Osvald
# ------------------------------------------------------------------

spectralnorm <- function(args) {
    n = if (length(args)) as.integer(args[[1]]) else 100L
    options(digits=10)

    eval_A <- function(i, j) 1 / ((i + j) * (i + j + 1) / 2 + i + 1)
    eval_A_times_u <- function(u) {
        ret <- double(n)
        for (i in 0:n1) {
            eval_A_col <- double(n)
            for (j in 0:n1)
	    eval_A_col[[j + 1]] <- eval_A(i, j)
            ret[[i + 1]] <- u %*% eval_A_col
        }
        return(ret)
    }
    eval_At_times_u <- function(u) {
        ret <- double(n)
        for (i in 0:n1) {
            eval_At_col <- double(n)
            for (j in 0:n1)
	    eval_At_col[[j + 1]] <- eval_A(j, i)
            ret[[i + 1]] <- u %*% eval_At_col
        }
        return(ret)
    }
    eval_AtA_times_u <- function(u) eval_At_times_u(eval_A_times_u(u))

    n1 <- n - 1
    u <- rep(1, n)
    v <- rep(0, n)
    for (itr in seq(10)) {
        v <- eval_AtA_times_u(u)
        u <- eval_AtA_times_u(v)
    }

    cat(sqrt(sum(u * v) / sum(v * v)), "\n")
}

if (!exists("i_am_wrapper"))
    spectralnorm(commandArgs(trailingOnly=TRUE))