clear all;
Sig = [2 1 0; 1 2 1; 0 1 2];

SigSqrt = [sqrt(2-sqrt(2))/2        -1  sqrt(2+sqrt(2))/2;
          -sqrt(2-sqrt(2))/sqrt(2)   0  sqrt(2+sqrt(2))/sqrt(2);
           sqrt(2-sqrt(2))/2         1  sqrt(2+sqrt(2))/2];
       
SigSqrtInv = ...
    [1/(2*sqrt(2-sqrt(2))) -1/(sqrt(2)*sqrt(2-sqrt(2))) 1/(2*sqrt(2-sqrt(2)));
    -1/2                    0                           1/2;
     1/(2*sqrt(2+sqrt(2)))  1/(sqrt(2)*sqrt(2+sqrt(2))) 1/(2*sqrt(2+sqrt(2)))];
 
% Part (b): Use Cholesky decomposition
SigSqrt = [sqrt(2)      0           0;
           1/sqrt(2)    sqrt(3/2)   0;
           0            sqrt(2/3)   sqrt(4/3)];
       
SigSqrtInv = [1/sqrt(2)      0           0;
             -1/sqrt(6)      sqrt(2/3)   0;
              1/(2*sqrt(3)) -1/sqrt(3)   sqrt(3)/2];