Skip to content
Prev 86757 / 398502 Next

Matrix indexing in a loop

There was an error in the function f -- it only worked correctly if
order was 1.  Here it is with that fixed.  The function f is the
only thing changed from my last post.  It makes use of the
fact that sum(x) == n and sum(x*x) == n*n  only occur when
one element of x is n and the rest are 0.

root3 <- function(x, order = 1) {
	f <- function(x,y) {
		z <- abs(as.numeric(x) - as.numeric(y))
		(sum(z) == order) & (sum(z*z) == order * order)
	}
	inner <- function(a,b,f)
			apply(b,1,function(x)apply(a,1,function(y)f(x,y)))

	Root <- function(x) {
		n <- length(dim(x))
		dd <- sapply(as.data.frame.table(x)[,-n-1], as.numeric)
		structure(inner(dd, dd, f) %*% c(x), .Dim = dim(x))
	}
	Root(x) / Root(1+0*x)
}

# tests
m <- matrix(1:24, 6) # 2d
root3(m)
mm <- array(1:24, 2:4)  # 3d
root3(mm)
On 2/17/06, Gabor Grothendieck <ggrothendieck at gmail.com> wrote: