🚀 go-pugleaf

   Some time ago I posted to this newsgroup, asking people to
email me matrix manipulation programs.  Here's what I received.
Thanks to everyone who responded (including those not detailed
below).
					  Simon

P.S.  Larry's programs will appear in the following book.
Lawrence C Paulson.
ML for the Working Programmer.
(Cambridge University Press, 1991, in press).
------------------------------------------------------------------






From jeroen%ruuinf.cs.ruu.nl@nsf.ac.uk Wed Dec  5 09:19:17 1990

This is my version, written in Miranda.
It stores matrices as lists of lists.
It avoids the list-indexing operator (!); instead
it uses functions like map, inits and tails as much as possible.

This satisfies my sense of "functional beauty", although
I must admit that functions like  inv  may be difficult to
read because of it.

--
Jeroen Fokker                                 |  jeroen@cs.ruu.nl
dept.of Computer Science, Utrecht University  |  tel.+31-30-534129
PObox 80089, 3508TB Utrecht, the Netherlands  |  fax.+31-30-513791
--



scalar == num
vector == [scalar]
matrix == [vector]

inprod :: vector -> vector -> scalar 
apply  :: matrix -> vector -> vector
mul    :: matrix -> matrix -> matrix
diag   :: matrix -> vector
transp :: matrix -> matrix
det    :: matrix -> scalar
inv    :: matrix -> matrix
solve  :: matrix -> vector -> vector

inprod a b = sum (map2 (*) a b)
apply  m v = map (inprod v) m

mul as bs = map f as
            where f a = map (inprod a)(transp bs)

diag []       =  []
diag (l:ls)   =  hd l :  diag (map tl ls)

antidiag = diag . reverse

transp [xs]     = map single xs
transp (xs:xss) = map2 (:) xs (transp xss)

gaps xs = map2 (++) (init(inits xs)) (tl(tails xs))
altern  = inprod plusmin

det [[x]]  =  x
det (v:m)  =  altern (map2 (*) v (map det (gaps (transp m))))

inv [[x]] = [[1/x]]
inv m = transp (zipp (*) c (mapp det (map gaps (map transp (gaps m)))))
        where c = mapp (/(det m)) (iterate tl plusmin)

solve m v = map ((/(det m)).det) (repls v (transp m))
            where repls v m = map2 (insert v) (init(inits m)) (tl(tails m))
                  insert v a b = a++ [v] ++b

m = [[1,3,1,1],[2,1,5,2],[1,-1,2,3],[4,1,-3,7]]
n=  [[2,1],[4,3]]

mapp f = map (map f)
zipp f = map2 (map2 f)

inits []     =  [[]]
inits (x:xs) =  [] : map (x:) (inits xs)

tails []     =  [[]]
tails (x:xs) =  (x:xs) : (tails xs)

plusmin = 1 : -1 : plusmin
single x = [x]

negvec = map neg
a $vecminus b = map2 (-) a b





From Larry.Paulson%uk.ac.cam.cl@cl.cam.ac.uk Wed Dec  5 10:37:03 1990

I enclose ML code for matrix transpose, multiplication, and Gaussian
Elimination.  It comes from my forthcoming book on Standard ML.  (The code is
entirely first-order because higher-order functions come later in the book!) 
Please report any errors you find; other comments are also welcome.  Hope you
find the code useful!
								Larry Paulson
------
(** Matrix transpose **)

fun headcol []    = []
  | headcol ((x::_) :: rows) = x :: headcol rows;

fun tailcols []    = []
  | tailcols ((_::xs) :: rows) = xs :: tailcols rows;

val matrix =
    [ ["a","b","c"],
      ["d","e","f"] ];

headcol matrix;
tailcols matrix;

fun transp ([]::rows) = []
  | transp rows = headcol rows :: transp (tailcols rows);

transp matrix;

(** Matrix product **)

fun dotprod([], []) = 0.0
  | dotprod(x::xs,y::ys) = x*y + dotprod(xs,ys);

fun rowprod(row, []) = []
  | rowprod(row, col::cols) =
        dotprod(row,col) :: rowprod(row,cols);

fun rowlistprod([], cols) = []
  | rowlistprod(row::rows, cols) =
        rowprod(row,cols) :: rowlistprod(rows,cols);

fun matprod(Arows,Brows) = rowlistprod(Arows, transp Brows);

val Arows = [ [2.0, 0.0], 
              [3.0, ~1.0],
              [0.0, 1.0],
              [1.0, 1.0] ];

val Brows = [ [1.0,  0.0, 2.0],
              [4.0, ~1.0, 0.0] ];

matprod(Arows,Brows);


(** Gaussian Elimination (Sedgewick Chapter 5) **)

(*the first row with absolutely greatest head*)
fun pivotrow [row] = row : real list
  | pivotrow (row1::row2::rows) =
      if abs(hd row1) >= abs(hd row2)
      then pivotrow(row1::rows)
      else pivotrow(row2::rows);

(*the matrix excluding the first row with head p*)
fun delrow (p, []) = []
  | delrow (p, row::rows) =
      if p = hd row then rows
      else row :: delrow(p, rows);

fun scalarprod(k, []) = [] : real list
  | scalarprod(k, x::xs) = k*x :: scalarprod(k,xs);

fun vectorsum ([], []) = [] : real list
  | vectorsum (x::xs,y::ys) = x+y :: vectorsum(xs,ys);

fun gausselim [row] = [row]
  | gausselim rows =
      let val p::prow = pivotrow rows
          fun elimcol [] = []
            | elimcol ((x::xs)::rows) =
                  vectorsum(xs, scalarprod(~x/p, prow))
                  :: elimcol rows
      in  (p::prow) :: gausselim(elimcol(delrow(p,rows)))
      end;

fun solutions [] = [~1.0]
  | solutions((x::xs)::rows) =
      let val solns = solutions rows
      in ~(dotprod(solns,xs)/x) :: solns  end;

(*Sedgewick's example*)
gausselim [[ 1.0,  3.0, ~4.0,  8.0], 
           [ 0.0,  2.0, ~2.0,  6.0], 
           [~1.0, ~2.0,  5.0, ~1.0]];
solutions it;
(*should be [1,5,2] with ~1 at end*)

gausselim [[ 0.0,  1.0,  2.0,  7.0,  7.0],
           [~4.0,  0.0,  0.0, ~5.0, ~17.0],
           [ 6.0, ~1.0,  0.0,  0.0,  28.0],
           [~2.0,  0.0,  2.0,  8.0,  12.0]];
solutions it;
(*should be [3,~10,5,1] with ~1 at end*)

------





From ahsnippe%praxis.cs.ruu.nl@nsf.ac.uk Fri Dec  7 08:20:10 1990

I've done some programming in Miranda which resulted in a library
for wire-frame-graphics. Part of it is a matrix-manipulation program.
A better version is owned by jeroen@cs.ruu.nl (Jeroen Fokker).


:
vector == [num]
matrix == [vector]



 dotprodukt [] [] = 0
 dotprodukt (a:l) (b:m) = a*b + dotprodukt l m

 apply xs y = map (dotprodukt y) xs

 ident n = map (line n) [0..(n-1)]   
            where line n i = (rep i 0)++[1]++(rep (n-i-1) 0)

 scalevec n x = map (*n) x

 scalemat n = map (scalevec n)

 mul x y = transpose (map (apply x) (transpose y))

 negvec = scalevec (-1)

 vecminus x y = map2 (-) x y

++++++++++++++++++++++++++

This is the one used in the graphics package

Next has some additonal functions :

++++++++++++++++++++++++++++++++

vector == [num]
matrix == [vector]

inprodukt :: vector -> vector -> num

 inprodukt [] [] = 0
 inprodukt (a:l) (b:m) = a*b + inprodukt l m


matapply :: matrix -> vector -> vector

 matapply xs y = map (inprodukt y) xs

transp :: matrix -> matrix

 transp m = (kop m) : (transp (rest m)) , if (vol m)
	 = []                          , otherwise

     vol []   = False
     vol (x:xs) = (x~=[])  
     
     kop  = map (hd) 
     
     rest = map (tl) 

ident :: num -> matrix

 ident n = map (line n) [0..(n-1)]   

 line n i = (rep i 0)++[1]++(rep (n-i-1) 0)


scalvec :: num -> vector -> vector

 scalvec n x = map (*n) x

scalmat :: num -> matrix -> matrix

 scalmat n = map (scalvec n)

matmul  :: matrix -> matrix -> matrix

 matmul x = map (matapply x)
det :: matrix -> num

 det  m = 0                      , if  ~(vol m) 
        = hd (hd m)              , if #m=1 
        = som (#m) m               , otherwise 

     som i m = (-1)^i * (((m)!0)!(i-1)) * (det (slicedmat i m)) + (som (i-1) m)  , if i >0
             = 0                                                             , otherwise

         slicedmat n []     = []
         slicedmat n (m:ms) = plak (stuk n ms) (rest2 n ms)
    
             stuk n []     = []
             stuk n (m:ms) = (take (n-1) m ) : stuk (n-1) ms 
            
             rest2 n []     = []
             rest2 n (m:ms) = (drop n m)  : rest2 n ms 
            
             plak [] [] = []
             plak (m:ms) (n:ns) = (m ++ n) : plak ms ns 

I hope the program is bug-free since I haven't used it for a long time.
If you have any questions, feel free to ask.


			  --[Arnold]--






From S.Clayman@cs.ucl.ac.uk Fri Dec  7 18:57:24 1990

I suspect you've had a bundle of replies. Anyway here's
a bit of Miranda I wrote ages ago
______________________________________________________________
|| some matrix things

%export scalar_mult_vector scalar_mult_matrix multmatrix

scalar == num
vector == [num]
matrix == [[num]]

|| this does s.v => [s.v1,s.v2 ...]
scalar_mult_vector::scalar->vector->vector
scalar_mult_vector s v = map ((*)s) v

|| this does s.M => [[s.m11,s.m12...],[s.m21,s.m22...]...]
scalar_mult_matrix::scalar->matrix->matrix
scalar_mult_matrix s m = map (scalar_mult_vector s) m

|| this is used in multmatrix to multiply
|| each row of matrix1 by matrix2
vector_mult_matrix::vector->matrix->vector
vector_mult_matrix v m
	= domult v m, #v=#m
	= error "vector_mult_matrix: vector and matrix have different lengths", otherwise
	  where
	  domult v ([]:any) = []
	  domult v m = sum (multi (map (*) v) (map hd m)) : domult v (map tl m)

multmatrix::matrix->matrix->matrix
multmatrix m1 m2 = multi (map vector_mult_matrix m1) (repeat m2)

|| multi - applies a list of functions to
|| a list of args 
|| multi [(*)1.5,(*)2,(+)3] [2,4,6]
|| should return [3,8,9]

multi::[*->**]->[*]->[**]
multi = map2 apply
	where
	apply f x = f x

______________________________________________________________

I hope it's not too obscure

Stuart






From pierre%ra.mcrcim.mcgill.edu@nsf.ac.uk Sat Dec  8 00:48:38 1990

I have done some work using MIT's eager functional language Id, to
perform LU decomposition.  My code is pretty simple-minded, and thus
inefficient, but code given by Arvind and Nikhil in a draft of a
forthcoming book on Id is far better written.  Perhaps they could be
convinced to mail you a copy.  One of their e-mail address is:

bug-id-world@mc.mit.edu

You can also look in J. Parallel and Distrib. Computing, vol 5 (no 5),
p. 460--493, for some of their earlier work.

		Pierre
--
------------------------------------------------------------------------------
 Simon Hughes, Dept. of Computing, Imperial College.  Email: sph@doc.ic.ac.uk