pfit <- function (x, y, z) {
  # empty object where to put returned quantities
  out <- NULL;     
  # design matrix
  out$X <- matrix( rep(1,length(x)), length(y), 1)
  out$X <- cbind(out$X, x, y)
  # inverse of (X^T)X, using Cholesky square matrix and
  # very efficient function crossprod(X) to calculate (X^T)X
  # [ the alternative way would have been
  #     inv.tX.X  <- solve( t(out$X) %*% out$X ]  
  inv.tX.X <- chol2inv( chol( crossprod(out$X) ))  
  # fit parameters
  out$par <- as.vector ( inv.tX.X %*% crossprod(out$X, z) )
  # residuals 
  out$res <- as.vector( z - out$X %*% out$par )
  # sum of squares of residuals
  out$res.srm <- sqrt(mean(out$res^2))
  # covariance matrix of fit parameters
  # including factor to take into account of the "n-n_p"
  out$cov <- inv.tX.X * out$res.srm^2 * length(x)/(length(x)- length(out$par) )
  # standard deviations
  out$std <- sqrt( diag( out$cov) )
  # correlation matrix
  out$rho <- out$cov / ( out$std %*% t(out$std) )
  return ( out )
}

dati <- read.table("dati_fit_piano.dat", header=TRUE)

# call hand-made function
r1 <- pfit(dati$x, dati$y, dati$z)

# call R function
r2 <- lm(dati$z ~ dati$x + dati$y)

# compare results
r1
summary(r2)