# chex35.r                         J F Monahan, August 2007
 # Example 3.5 Condition of Hilbert matrix
 #
 # set n = size of matrix
  n <- 5                          # a nice size
  "size of matrix"
  n
 # Hilbert matrix is easy, but not exactly represented
  A <- matrix(0,n,n)              # initialize to zero
    for (i in 1:n) {              # fill element by element
      for (j in 1:n) { A[i,j] <- 1/(i+j-1) } }
  "Hilbert matrix"
  A
 # inverse of matrix has integer values -- exact representation
 # save some constants to make calculation easier
  Is <- 1:n
  bi <- choose(Is+n-1,Is-1)       # save these (use with i)
  bj <- choose(n,Is)              # and these (use with j)
  B <- matrix(0,n,n)              # initialize to zero
    sign <- -1
    for (i in 1:n) {
      for (j in 1:n) {
        sign <- -sign              # (-1)**(i+j)
        B[i,j] <- sign*j*choose(i+j-2,i-1)*bi[i]*choose(j+n-1,n-i)*bj[j]
                      } }   
  "inverse of Hilbert matrix"
  B
 # is it really the inverse
  "A %*% B"
  A %*% B
 # compute p = inf norms (row sum norm)
  normA <- max( apply(abs(A),1,sum) )
  "infinity norm of Hilbert matrix"
  normA
  normB <- max( apply(abs(B),1,sum) )
  "infinity norm of inverse of Hilbert matrix"
  normB
  "kappa and reciprocal"
  kappa <- normA*normB
  kappa
  1/kappa
 # 
  rm(list=ls())