function(X, k = 0, sca = "A", scores = F)
{
#
# Robust Principal Component Analysis
#
# C. Croux (Brussels, Belgium)
#
# k is the  number of components we want
# sca is a scale estimator
#
	n <- nrow(X)
	p <- ncol(X)
	if(k == 0)
		p1 <- min(n, p)
	else p1 <- k
	S <- rep(1, p1)
	V <- matrix(1:(p * p1), ncol = p1, nrow = p)
	P <- diag(p)	#m <- l1median(X)	
	m <- apply(X, 2, median)
	Xcentr <- scale(X, center = m, scale = F)
	for(k in 1:p1) {
		B <- Xcentr %*% P
		Bnorm <- apply(B, 1, vecnorm)
		A <- diag(1/Bnorm) %*% B
		Y <- A %*% P %*% t(X)
		if(sca == "mad")
			s <- apply(Y, 1, mad)
		if(sca == "tau")
			s <- apply(Y, 1, scale.tau)
		if(sca == "A")
			s <- apply(Y, 1, scale.a)
		j <- order(s)[n]
		S[k] <- s[j]
		V[, k] <- A[j,  ]
		if(V[1, k] < 0)
			V[, k] <- (-1) * V[, k]
		P <- P - (V[, k] %*% t(V[, k]))
	}
	if(scores) {
		list(scale = S, rotation = V, scores = Xcentr %*% V)
	}
	else list(scale = S, rotation = V)
}