R : Copyright 2005, The R Foundation for Statistical Computing
Version 2.1.1  (2005-06-20), ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for a HTML browser interface to help.
Type 'q()' to quit R.

>  # vadjust.r                        J F Monahan, September 2008
>  #   example of vadjust -- vector version of adjust
>  #
>  #                  simple matrix -- scaled second column
>  X <- matrix(c(rep(1,6),(1:6)/10,(1:6)^2),6,3)
>  "design matrix"
[1] "design matrix"
>  X
     [,1] [,2] [,3]
[1,]    1  0.1    1
[2,]    1  0.2    4
[3,]    1  0.3    9
[4,]    1  0.4   16
[5,]    1  0.5   25
[6,]    1  0.6   36
>  "inner product -- simple matrix"
[1] "inner product -- simple matrix"
>  xpx <- t(X) %*% X
>  xpx
     [,1]  [,2]   [,3]
[1,]  6.0  2.10   91.0
[2,]  2.1  0.91   44.1
[3,] 91.0 44.10 2275.0
>  Lt <- chol(xpx)
>  "Cholesky factor -- transposed to be upper triangular"
[1] "Cholesky factor -- transposed to be upper triangular"
>  Lt
         [,1]      [,2]      [,3]
[1,] 2.449490 0.8573214 37.150594
[2,] 0.000000 0.4183300 29.283101
[3,] 0.000000 0.0000000  6.110101
>  det <-prod(diag(Lt))
>  "determinant"
[1] "determinant"
>  det*det
[1] 39.2
>  # use function 'adjust' previously constructed
>   adjust <- function(x) {
+      if( x <= 0 ) c( 0, -2**31 ) else {  # flag if not positive
+           i <- 0
+           while (x > 16 ) { 
+              x <- x/16
+              i <- i + 4
+                           }
+           while (x < 1 ) { 
+              x <- x*16
+              i <- i - 4
+                           }
+           c(x,i)
+                                       }
+                                       }
>  # now compute -- it may have already overflowed
>  adjust(det*det)
[1] 2.45 4.00
>  #
>  # modification of function 'adjust' previously constructed
>  #  now vector version
>   vadjust <- function(x) {
+      i <- rep(0,length(x))
+      zneg <- ( x <= 0 )       # flag if not positive
+      i[zneg] <- -(2^31)
+      x[zneg] <- 0
+      big <- (x > 2)           # different scaling, 16 may be better
+           while ( max(big) > 0 ) { 
+              x[big] <- x[big]/2         # divide
+              i[big] <- i[big] + 1       # adjust exponent
+              big <- (x > 2)
+                           }
+      small <- (x < 1)
+           while ( max(small) > 0 ) { 
+              x[small] <- x[small]*2     # multiply
+              i[small] <- i[small] - 1   # adjust exponent
+              small <- (x < 1)
+                           }
+           list(x,i)       }                 # make list of vectors
>  #
>  # do again in a more stable fashion
>  both <- vadjust(diag(Lt))  
>  "both is a list"
[1] "both is a list"
>  both
[[1]]
[1] 1.224745 1.673320 1.527525

[[2]]
[1]  1 -2  2

>  prod(both[[1]])              # part of the determinant
[1] 3.130495
> 
>  sum(both[[2]])               # exponent part
[1] 1
>  "determinant"
[1] "determinant"
>  ((prod(both[[1]]))^2)*(2^(2*sum(both[[2]])))
[1] 39.2
>  #                            done
>  rm(list=ls())
>  q()