Algorithm Implementation/Linear Algebra/Determinant of a Matrix

From Wikibooks, open books for an open world
Jump to: navigation, search

Haskell[edit]

--Takes input with the matrix as a list of rows, throws an error if the list isn't square
--This probably isn't the most elegant/efficient way to do it, but it's something
alternatingSum :: Num a => [a] -> a
alternatingSum [] = 0
alternatingSum l  = iter l 0 0
    where iter (x:xs) n result = if even n then iter xs (n + 1) (result + x) else iter xs (n + 1) (result - x)
          iter [] _ result = result
 
removeColumn :: Int -> [[a]] -> [[a]]
removeColumn n = map (removeIndex n)
 
removeIndex :: Int -> [a] -> [a]
removeIndex 0 (x:xs) = xs
removeIndex n (x:xs) = x : removeIndex (n-1) xs
removeIndex _ []     = []
 
isSquare :: [[a]] -> Bool
isSquare (row:rows) = and $ map eqFirstLength rows
    where l = length row
          eqFirstLength list = length list == l
isSquare [] = True
 
determinant :: Num a => [[a]] -> a
determinant matrix 
    | isSquare matrix   = helper matrix
    | otherwise         = error "DETERMINANT --> MATRIX IS NOT SQUARE"
    where helper [[x1,x2],[y1,y2]] = x1 * y2 - x2 * y1
          helper (x:xs)            = let l = length x in alternatingSum $ zipWith (*) x $ map (\n -> determinant $ removeColumn n xs) [0..l]
          helper [] = 0

Java[edit]

       public static double[][] reduce(double[][] x, double[][] y, int r, int c, int n)
      {
         for (int h = 0, j = 0; h < n; h++)
         {
            if (h == r) 
               continue;
            for (int i = 0, k = 0; i < n; i++)
            {
               if (i == c) 
                  continue;
               y[j][k] = x[h][i];
               k++;
            } //end inner loop
            j++;
         } //end outer loop
         return y;
      } //end method
   //===================================================
       public static double det(int NMAX, double[][] x)
      {
         double ret=0;
         if (NMAX < 4)//base case
         {
            double prod1=1, prod2=1;
            for (int i = 0; i < NMAX; i++)
            {
                  prod1=1;
                  prod2=1;
 
               for (int j = 0; j < NMAX; j++)
               {
                  prod1 *= x[(j + i + 1) % NMAX][j];
                  prod2 *= x[(j + i + 1) % NMAX][NMAX - j - 1];
               } //end inner loop
               ret += prod1 - prod2;
            } //end outer loop
            return ret;
         } //end base case
         double[][] y = new double [NMAX - 1][NMAX - 1];
         for (int h = 0; h < NMAX; h++)
         {
            if (x[h][0] == 0) 
               continue;
            reduce(x, y, h, 0, NMAX);
            if (h % 2 == 0) ret -= det(NMAX - 1, y) * x[h][0];
            if (h  %2 == 1) ret += det(NMAX - 1, y) * x[h][0];
         } //end loop
         return ret;
      } //end method