Cramer's rule: Difference between revisions

Content added Content deleted
Line 419: Line 419:
Just [2.0000000000000004,-11.999999999999998,-4.0,0.9999999999999999]
Just [2.0000000000000004,-11.999999999999998,-4.0,0.9999999999999999]
===Version 2===
===Version 2===
We use Rational numbers for more precision.
a % b is the rational a / b.
<lang Haskell>
<lang Haskell>
s_permutations :: [a] -> [([a], Int)]
s_permutations :: [a] -> [([a], Int)]
Line 427: Line 429:
insertEv x [] = [[x]]
insertEv x [] = [[x]]
insertEv x l@(y:ys) = (x:l) : map (y:) (insertEv x ys)
insertEv x l@(y:ys) = (x:l) : map (y:) (insertEv x ys)
foldlZipWith::(a -> b -> c) -> (d -> c -> d) -> d -> [a] -> [b] -> d
foldlZipWith::(a -> b -> c) -> (d -> c -> d) -> d -> [a] -> [b] -> d
foldlZipWith _ _ u [] _ = u
foldlZipWith _ _ u [] _ = u
foldlZipWith _ _ u _ [] = u
foldlZipWith _ _ u _ [] = u
foldlZipWith f g u (x:xs) (y:ys) = foldlZipWith f g (g u (f x y)) xs ys
foldlZipWith f g u (x:xs) (y:ys) = foldlZipWith f g (g u (f x y)) xs ys

foldl1ZipWith::(a -> b -> c) -> (c -> c -> c) -> [a] -> [b] -> c
foldl1ZipWith::(a -> b -> c) -> (c -> c -> c) -> [a] -> [b] -> c
foldl1ZipWith _ _ [] _ = error "First list is empty"
foldl1ZipWith _ _ [] _ = error "First list is empty"
foldl1ZipWith _ _ _ [] = error "Second list is empty"
foldl1ZipWith _ _ _ [] = error "Second list is empty"
foldl1ZipWith f g (x:xs) (y:ys) = foldlZipWith f g (f x y) xs ys
foldl1ZipWith f g (x:xs) (y:ys) = foldlZipWith f g (f x y) xs ys

multAdd::(a -> b -> c) -> (c -> c -> c) -> [[a]] -> [[b]] -> [[c]]
multAdd::(a -> b -> c) -> (c -> c -> c) -> [[a]] -> [[b]] -> [[c]]
multAdd f g xs ys = map (\us -> foldl1ZipWith (\u vs -> map (f u) vs) (zipWith g) us ys) xs
multAdd f g xs ys = map (\us -> foldl1ZipWith (\u vs -> map (f u) vs) (zipWith g) us ys) xs

mult:: Num a => [[a]] -> [[a]] -> [[a]]
mult xs ys = multAdd (*) (+) xs ys
elemPos::[[a]] -> Int -> Int -> a
elemPos ms i j = (ms !! i) !! j
prod:: Num a => ([[a]] -> Int -> Int -> a) -> [[a]] -> [Int] -> a
prod f ms = product.zipWith (f ms) [0..]
s_determinant:: Num a => ([[a]] -> Int -> Int -> a) -> [[a]] -> [([Int],Int)] -> a
s_determinant f ms = sum.map (\(is,s) -> fromIntegral s * prod f ms is)
determinant:: Num a => [[a]] -> a
determinant ms = s_determinant elemPos ms $ s_permutations [0..pred.length $ ms]
elemCramerPos::Int -> [[a]] -> [[a]] -> Int -> Int -> a
elemCramerPos l ks ms i j = if j /= l then elemPos ms i j else head (ks !! i)
solveCramer::(Fractional a, Eq a) => [[a]] -> [[a]] -> [[a]]
solveCramer ms ks = xs
where
xs | d /= 0 = map (\l -> return.(/d) $ s_determinant (elemCramerPos l ks) ms ps) $ ls
| otherwise = []
ls = [0..pred.length $ ms]
ps = s_permutations ls
d = s_determinant elemPos ms ps
task = do
let a = [[2,-1, 5, 1]
,[3, 2, 2,-6]
,[1, 3, 3,-1]
,[5,-2,-3, 3]]
let y = [[-3], [-32], [-47], [49]]
let x = solveCramer a y
let b = mult a x
putStrLn "a ="
mapM_ print a
putStrLn "y ="
mapM_ print y
putStrLn "solve: a * x = y => x = solveCramer a y ="
mapM_ print x
putStrLn "verification: b = a * x ="
mapM_ print b
putStrLn $ "test: b == y = "
print $ b == y
main = task
</lang>
{{out}}
<pre>
a =
[2.0,-1.0,5.0,1.0]
[3.0,2.0,2.0,-6.0]
[1.0,3.0,3.0,-1.0]
[5.0,-2.0,-3.0,3.0]
y =
[-3.0]
[-32.0]
[-47.0]
[49.0]
solve: a * x = y => x = solveCramer a y =
[2.0]
[-12.0]
[-4.0]
[1.0]
verification: b = a * x =
[-3.0]
[-32.0]
[-47.0]
[49.0]
test: b == y =
True
</pre>
===Version 3===
Almost the same than Version 2 but more general and we use Rational numbers for more precision.
a % b is the rational a / b.
<lang Haskell>
s_permutations :: [a] -> [([a], Int)]
s_permutations = flip zip (cycle [1, -1]) . (foldl aux [[]])
where aux items x = do
(f,item) <- zip (cycle [reverse,id]) items
f (insertEv x item)
insertEv x [] = [[x]]
insertEv x l@(y:ys) = (x:l) : map (y:) (insertEv x ys)
foldlZipWith::(a -> b -> c) -> (d -> c -> d) -> d -> [a] -> [b] -> d
foldlZipWith _ _ u [] _ = u
foldlZipWith _ _ u _ [] = u
foldlZipWith f g u (x:xs) (y:ys) = foldlZipWith f g (g u (f x y)) xs ys

foldl1ZipWith::(a -> b -> c) -> (c -> c -> c) -> [a] -> [b] -> c
foldl1ZipWith _ _ [] _ = error "First list is empty"
foldl1ZipWith _ _ _ [] = error "Second list is empty"
foldl1ZipWith f g (x:xs) (y:ys) = foldlZipWith f g (f x y) xs ys

multAdd::(a -> b -> c) -> (c -> c -> c) -> [[a]] -> [[b]] -> [[c]]
multAdd f g xs ys = map (\us -> foldl1ZipWith (\u vs -> map (f u) vs) (zipWith g) us ys) xs

mult:: Num a => [[a]] -> [[a]] -> [[a]]
mult:: Num a => [[a]] -> [[a]] -> [[a]]
mult xs ys = multAdd (*) (+) xs ys
mult xs ys = multAdd (*) (+) xs ys
Line 550: Line 454:
s_determinant:: Num a => ([[a]] -> Int -> Int -> a) -> [[a]] -> [([Int],Int)] -> a
s_determinant:: Num a => ([[a]] -> Int -> Int -> a) -> [[a]] -> [([Int],Int)] -> a
s_determinant f ms = sum.map (\(is,s) -> fromIntegral s * prod f ms is)
s_determinant f ms = sum.map (\(is,s) -> fromIntegral s * prod f ms is)
determinant:: Num a => [[a]] -> a
determinant ms = s_determinant elemPos ms $ s_permutations [0..pred.length $ ms]
elemCramerPos::Int -> Int -> [[a]] -> [[a]] -> Int -> Int -> a
elemCramerPos::Int -> Int -> [[a]] -> [[a]] -> Int -> Int -> a
Line 568: Line 469:
ps = s_permutations ls
ps = s_permutations ls
d = s_determinant elemPos ms ps
d = s_determinant elemPos ms ps
matI::(Num a) => Int -> [[a]]

matI n = [ [fromIntegral.fromEnum $ i == j | i <- [1..n]] | j <- [1..n]]
matI n = [ [fromIntegral.fromEnum $ i == j | i <- [1..n]] | j <- [1..n]]

task::[[Rational]] -> [[Rational]] -> IO()
task = do
let a = [[2,-1, 5, 1]
task a b = do
,[3, 2, 2,-6]
let x = solveCramer a b
,[1, 3, 3,-1]
let u = map (map fromRational) x
,[5,-2,-3, 3]]
let y = mult a x
let y = [[-3], [-32], [-47], [49]]
let identity = matI (length x)
let x = solveCramer a y
let a1 = solveCramer a identity
let b = mult a x
let h = mult a a1
let i4 = matI 4
let z = mult a1 b
let a1 = solveCramer a i4
let h4 = mult a a1
let z = mult a1 y
putStrLn "a ="
putStrLn "a ="
mapM_ print a
mapM_ print a
putStrLn "y ="
putStrLn "b ="
mapM_ print y
mapM_ print b
putStrLn "solve: a * x = y => x = solveCramer a y ="
putStrLn "solve: a * x = b => x = solveCramer a b ="
mapM_ print x
mapM_ print x
putStrLn "verification: b = a * x ="
putStrLn "u = fromRationaltoDouble x ="
mapM_ print b
mapM_ print u
putStrLn $ "test: b == y = "
putStrLn "verification: y = a * x = mult a x ="
print $ b == y
mapM_ print y
putStrLn "identity matrix: i4 ="
putStrLn $ "test: y == b = "
mapM_ print i4
print $ y == b
putStrLn "find: a1 = inv(a) => solve: a * a1 = i4 => a1 = solveCramer a i4 ="
putStrLn "identity matrix: identity ="
mapM_ print identity
putStrLn "find: a1 = inv(a) => solve: a * a1 = identity => a1 = solveCramer a identity ="
mapM_ print a1
mapM_ print a1
putStrLn "verification: h4 = a * a1 ="
putStrLn "verification: h = a * a1 = mult a a1 ="
mapM_ print h4
mapM_ print h
putStrLn $ "test: h4 == i4 = "
putStrLn $ "test: h == identity = "
print $ h4 == i4
print $ h == identity
putStrLn "z = a1 * y ="
putStrLn "z = a1 * b = mult a1 b ="
mapM_ print z
mapM_ print z
putStrLn "test: z == x ="
putStrLn "test: z == x ="
print $ z == x
print $ z == x
main = task
main = do
let a = [[2,-1, 5, 1]
,[3, 2, 2,-6]
,[1, 3, 3,-1]
,[5,-2,-3, 3]]
let b = [[-3], [-32], [-47], [49]]
task a b
</lang>
</lang>
{{out}}
{{out}}
Line 615: Line 521:
[1 % 1,3 % 1,3 % 1,(-1) % 1]
[1 % 1,3 % 1,3 % 1,(-1) % 1]
[5 % 1,(-2) % 1,(-3) % 1,3 % 1]
[5 % 1,(-2) % 1,(-3) % 1,3 % 1]
y =
b =
[(-3) % 1]
[(-3) % 1]
[(-32) % 1]
[(-32) % 1]
[(-47) % 1]
[(-47) % 1]
[49 % 1]
[49 % 1]
solve: a * x = y => x = solveCramer a y =
solve: a * x = b => x = solveCramer a b =
[2 % 1]
[2 % 1]
[(-12) % 1]
[(-12) % 1]
[(-4) % 1]
[(-4) % 1]
[1 % 1]
[1 % 1]
u = fromRationaltoDouble x =
verification: b = a * x =
[2.0]
[-12.0]
[-4.0]
[1.0]
verification: y = a * x = mult a x =
[(-3) % 1]
[(-3) % 1]
[(-32) % 1]
[(-32) % 1]
[(-47) % 1]
[(-47) % 1]
[49 % 1]
[49 % 1]
test: b == y =
test: y == b =
True
True
identity matrix: i4 =
identity matrix: identity =
[1 % 1,0 % 1,0 % 1,0 % 1]
[1 % 1,0 % 1,0 % 1,0 % 1]
[0 % 1,1 % 1,0 % 1,0 % 1]
[0 % 1,1 % 1,0 % 1,0 % 1]
[0 % 1,0 % 1,1 % 1,0 % 1]
[0 % 1,0 % 1,1 % 1,0 % 1]
[0 % 1,0 % 1,0 % 1,1 % 1]
[0 % 1,0 % 1,0 % 1,1 % 1]
find: a1 = inv(a) => solve: a * a1 = i4 => a1 = solveCramer a i4 =
find: a1 = inv(a) => solve: a * a1 = identity => a1 = solveCramer a identity =
[4 % 171,11 % 171,10 % 171,8 % 57]
[4 % 171,11 % 171,10 % 171,8 % 57]
[(-55) % 342,(-23) % 342,119 % 342,2 % 57]
[(-55) % 342,(-23) % 342,119 % 342,2 % 57]
[107 % 684,(-5) % 684,11 % 684,(-7) % 114]
[107 % 684,(-5) % 684,11 % 684,(-7) % 114]
[7 % 684,(-109) % 684,103 % 684,7 % 114]
[7 % 684,(-109) % 684,103 % 684,7 % 114]
verification: h4 = a * a1 =
verification: h = a * a1 = mult a a1 =
[1 % 1,0 % 1,0 % 1,0 % 1]
[1 % 1,0 % 1,0 % 1,0 % 1]
[0 % 1,1 % 1,0 % 1,0 % 1]
[0 % 1,1 % 1,0 % 1,0 % 1]
[0 % 1,0 % 1,1 % 1,0 % 1]
[0 % 1,0 % 1,1 % 1,0 % 1]
[0 % 1,0 % 1,0 % 1,1 % 1]
[0 % 1,0 % 1,0 % 1,1 % 1]
test: h4 == i4 =
test: h == identity =
True
True
z = a1 * y =
z = a1 * b = mult a1 b =
[2 % 1]
[2 % 1]
[(-12) % 1]
[(-12) % 1]
Retrieved from "https://rosettacode.org/wiki/Cramer%27s_rule"