flatten [] = [] flatten (x : xs) = x ++ (flatten xs) -- ### I/O ### ---------------------------------------------------------------- -- reads first n lines from stdin and returns them in a list getLines :: (Ord t, Num t) => t -> IO [String] getLines n | n <= 0 = return [] | otherwise = do x <- getLine xs <- getLines (n - 1) return (x:xs) -- transforms a list of lists to a printable matrix string vec_to_str [] = "\n" vec_to_str (x : xs) = (take 6 ((if x < 0 then "" else " ") ++ (show x) ++ [' ' | _ <- [1..]])) ++ "\t" ++ (vec_to_str xs) mat_to_str [] = "" mat_to_str (row : rows) = (vec_to_str row) ++ (mat_to_str rows) -- given a string representing vector b, finds vector x s.t. Ax = b get_solution_string str_b lower upper swaps = let b = map read (words str_b) :: [Double] y = forward_sub 1 lower [] (apply_swaps b swaps) x = backward_sub upper y in "Ax = " ++ (show b) ++ " => x = " ++ (show x) ++ "\n" -- given a list of strings of vectors, returns a string of solutions get_all_solutions strings lower upper swaps = flatten [get_solution_string b lower upper swaps | b <- strings] -- ######################################################################### -- -- ### Linear Algebra ### ----------------------------------------------------- mul_vv [] [] = 0 mul_vv (x : xs) (y : ys) = x * y + (mul_vv xs ys) div_vc [] _ = [] div_vc (x : xs) c = (x / c) : (div_vc xs c) div_mc [] _ = [] div_mc (row : rows) c = (div_vc row c) : (div_mc rows c) outer_product [] _ = [] outer_product (x : xs) ys = [x * y | y <- ys] : (outer_product xs ys) sub_vv [] [] = [] sub_vv (x : xs) (y : ys) = (x - y) : (sub_vv xs ys) sub_mm [] [] = [] sub_mm (row1 : mat1) (row2 : mat2) = (sub_vv row1 row2) : (sub_mm mat1 mat2) -- ######################################################################### -- -- ### Mics ### ---------------------------------------------------------------- -- splits the matrix returning the first column and the rest of the matrix split [] = ([], []) split ((x : xs) : rows) = let (col, mat) = split rows in ((x : col), (xs : mat)) -- reverse of split merge [] [] = [] merge (x : xs) (row : mat) = ((x : row) : (merge xs mat)) -- swaps two elements of list at positions i1, j1 swap i1 j1 list = if i1 == j1 then list else let i = minimum [i1, j1] j = maximum [i1, j1] first = take i list second = drop (i + 1) (take j list) third = drop (j + 1) list in first ++ [list !! j] ++ second ++ [list !! i] ++ third apply_swaps list [] = list apply_swaps list ((i, j) : swaps) = apply_swaps (swap i j list) swaps -- ######################################################################### -- -- ### LUP decomposition ### --------------------------------------------------- -- returns the index of the maximum element in the list max_idx _ [] _ res = res max_idx i (x : xs) cur_max res = if (abs x) > cur_max then max_idx (i + 1) xs (abs x) i else max_idx (i + 1) xs cur_max res -- pivots the matrix bringing the biggest element to the top left corner. -- returns the new matrix and the permutation represented as a lsit of swaps. pivot mat = let (row, rest) = split mat i = max_idx 0 row (-1) (-1) in ((swap 0 i mat), (0, i)) add_1 [] = [] add_1 ((i, j) : rest) = (((i + 1), (j + 1)) : (add_1 rest)) -- computes the LUP decomposition returning (L, U, P). lup_dec [[x]] = ([[1]], [[x]], []) lup_dec input_mat = let (((a : row) : rows), q) = pivot input_mat (col, mat) = split rows schur = sub_mm mat (div_mc (outer_product col row) a) (l1, u1, p1) = lup_dec schur (swap_i, swap_j) = head p1 l_top = (1 : [0 | _ <- row]) u_left = (a : [0 | _ <- col]) l = (l_top : (merge (div_vc (apply_swaps col p1) a) l1)) u = merge u_left (row : u1) p = (q : (add_1 p1)) in (l, u, p) --in ((a : row) : rows, schur, [q]) -- ######################################################################### -- -- ### Determinant ### -------------------------------------------------------- determinant [] = 1 determinant ((x : row) : upper) = let (col, mat) = split upper in x * (determinant mat) -- ######################################################################### -- -- ### Linear equations ### --------------------------------------------------- forward_sub _ [] xs _ = xs forward_sub i (row : lower) xs (b : bs) = let prefix = mul_vv [x | (x, j) <- (zip row [1..(i - 1)])] xs x = b - prefix in forward_sub (i + 1) lower (xs ++ [x]) bs backward_sub [[x]] [b] = [b / x] backward_sub ((x : row) : upper) (b : bs) = let (col, mat) = split upper xs = backward_sub mat bs prefix = mul_vv row xs in ((b - prefix) / x) : xs -- ######################################################################### -- main :: IO() main = do -- reading and parsing the size of the matrix n_s <- getLine let n = read n_s :: Integer -- reading rows of the matrix as a list of strings rows <- getLines n -- transforming the list of strings to a double array of Integers let mat = [row | row_s <- rows, let row = (map read (words row_s) :: [Double])] -- calculating the LU-decomposition let (lower, upper, p) = lup_dec mat ---- calculating the determinant using the decomposed matrices let det = determinant upper * (if ((length p) `div` 2) == 0 then 1 else -1) putStrLn $ "Input matrix: \n" ++ (mat_to_str mat) ++ "\n\n" putStrLn $ "Unit-Lower Triangular: \n" ++ (mat_to_str lower) ++ "\n" putStrLn $ "Upper Triangular: \n" ++ (mat_to_str upper) ++ "\n" putStrLn $ "Permutations: " ++ (show p) ++ "\n" putStrLn $ "Determinant: " ++ (show det) ++ "\n" -- reading the amount of queries for linear equation solutions k_s <- getLine let k = read k_s :: Integer -- reading queries as a list of strings queries <- getLines k let solutions = get_all_solutions queries lower upper p putStrLn $ "\n\nSolutions: \n\n" ++ solutions