Many times in this list people said that Haskell is in-efficient in
comparison to non-functional tools, like C and others.
And often refer to the matrix operations, linear algebra tasks,
saying how the arrays are useful and how they need to be mutable.
Maybe we arrange some simple and clear benchmark?
The enclosed programs are written boldly in Haskell for the two ways
to compute the determinant.
A matrix is represented as a list of lists. No arrays. Nothing is
mutable "in place", God save.
To avoid large numbers, we put the number domain C to be integers
modulo prime p.
cNeg, cAdd ... cDiv program explicitly the arithmetic modulo p.
Now, could someone, please, write similar programs in other
languages, say, C ?
Use mutable arrays, as you like, and so on.
Present the source program, please.
Then, we will run det m, det_G m
for the matrix of size n = 20, 40 ... and see the timing.
?
------------------
Sergey Mechveliani
[EMAIL PROTECTED]
--------------------------------------------------------------------
det_G :: [[C]] -> C
-- determinant of square matrix computed by Gauss method
det_G [[a]] = a
det_G rows = case span ((==0).head) rows of
(_ , [] ) -> 0 -- zero column
(zhs, (a:row):rows') ->
(case
map (\ (b:row') -> mulSubRow (cDiv b a) row' row) rows'
of
rs -> makeSign zhs $ cMul a $ det_G ((map tail zhs)++rs)
)
where
makeSign (_:_:xs) = makeSign xs
makeSign [] = id
makeSign [_] = cNeg
mulSubRow q = zipWith (\b a -> cSub b $ cMul q a)
--------------------------------------------------------------------
det :: [[C]] -> C -- determinant of square matrix computed by
-- expanding by the row
det [[a]] = a
det (row:rows) = dt (1::Int) '+' row rows (0 :: C)
-- initiate det = 0, sign= '+',
-- for i = 1 to n ...
where
dt _ _ [] _ res = res
dt i sign (a:row) rows res =
if
a==0 then dt (i+1) (invSign sign) row rows res --skip zero
else
(case (sign, det $ delColumn i rows)
of
('+',minorDet) -> recurse '-' cAdd minorDet
('-',minorDet) -> recurse '+' cSub minorDet
)
where
recurse s addOrSub minorDet =
dt (i+1) s row rows (addOrSub res (cMul a minorDet))
invSign '+' = '-'
invSign '-' = '+'
delColumn :: Int -> [[a]] -> [[a]]
delColumn n = map (del_n_th n)
del_n_th :: Int -> [a] -> [a] -- remove element No n from list
del_n_th _ [] = []
del_n_th 0 xs = xs
del_n_th 1 (_:xs) = xs
del_n_th n (x:xs) = x:(del_n_th (n-1) xs)
----------------------------------------------------------------------
type C = Int -- base domain - only it is viewed modulo p
p = 3 :: C
cNeg 0 = 0
cNeg a = p-a :: C
cAdd, cSub, cMul, cDiv :: C -> C -> C
cAdd a b = mod (a+b) p
cSub a b = cAdd a $ cNeg b
cMul a b = mod (a*b) p
cPow :: C -> Int -> C
cPow a n = mod (a^n) p
cDiv a b = cMul a $ cInv b
cInv 0 = error "cInv 0"
cInv a = case gcdE a p
of
(1,u,v) -> case mod u p
of
u' -> if u' < 0 then u'+p else u'
_ -> error "gcdE a p should yeild (1,u,v)"
gcdE :: C -> C -> (C,C,C) --extended GCD: g = gcd(x,y) = u*x+v*y
--x y g u v
gcdE x y = gc (1,0,x) (0,1,y)
-- (u1,u2,u3) (v1,v2,v3)
-- Euclidean GCD is applied to u3,v3, operations on
-- u1,u2, v1,v2 are parallel to ones of u3,v3
where
gc (u1,u2,u3) (v1,v2,v3) =
if v3==0 then (u3,u1,u2)
else
case quotRem u3 v3
of
(q,r) -> gc (v1,v2,v3) (u1-q*v1, u2-q*v2, r)
