open System
open Microsoft.FSharp.Collections

type Matrix = class
   val m : float[,]
   new (rows : int, cols : int) =
       {
           m = Array2D.create rows cols 0.0
       }
   member x.Item
       with get (i, j) = x.m.[i - 1, j - 1]
       and set (i, j) v = x.m.[i - 1, j - 1] <- v
   member x.Get (i : int) (j : int) = x.m.[i - 1, j - 1]
   member x.Set (i : int) (j : int) (v : float) = x.m.[i - 1, j - 1] <- v
   member x.Rows
       with get () = x.m.GetLength(0)
   member x.Cols
       with get () = x.m.GetLength(1)
   member x.T () : Matrix =
       let ret = new Matrix(x.Cols, x.Rows)
       for i = 0 to x.Rows - 1 do
           for j = 0 to x.Cols - 1 do
               ret.m.[j, i] <- x.m.[i, j]
       ret
   static member Identity (n : int) : Matrix =
       let ret = new Matrix(n, n)
       for i = 0 to n - 1 do
           ret.m.[i, i] <- 1.0
       ret
   static member (+) (x : Matrix, y : Matrix) : Matrix =
       if x.Rows <> y.Rows || x.Cols <> y.Cols then raise (new Exception())
       let ret = new Matrix(x.Rows, x.Cols)
       for i = 0 to ret.Rows - 1 do
           for j = 0 to ret.Cols - 1 do
               ret.m.[i, j] <- x.m.[i, j] + y.m.[i, j]
       ret
   static member (-) (x : Matrix, y : Matrix) : Matrix =
       if x.Rows <> y.Rows || x.Cols <> y.Cols then raise (new Exception())
       let ret = new Matrix(x.Rows, x.Cols)
       for i = 0 to ret.Rows - 1 do
           for j = 0 to ret.Cols - 1 do
               ret.m.[i, j] <- x.m.[i, j] - y.m.[i, j]
       ret
   static member (*) (k : float, x : Matrix) : Matrix =
       let ret = new Matrix(x.Rows, x.Cols)
       for i = 0 to ret.Rows - 1 do
           for j = 0 to ret.Cols - 1 do
               ret.m.[i, j] <- k * x.m.[i, j]
       ret
   static member (*) (x : Matrix, k : float) : Matrix = k * x
   static member (/) (x : Matrix, k : float) : Matrix = (1.0 / k) * x
   static member (*) (x : Matrix, y : Matrix) : Matrix =
       if x.Cols <> y.Rows then raise (new Exception())
       let M = x.Rows
       let S = x.Cols
       let N = y.Cols
       let ret = new Matrix(M, N)
       for i = 0 to M - 1 do
           for j = 0 to N - 1 do
               for k = 0 to S - 1 do
                   ret.m.[i, j] <- ret.m.[i, j] + x.m.[i, k] * y.m.[k, j]
       ret
   static member Dot (x : Matrix) (y : Matrix) : float =
       if (x.Rows <> 1 && x.Cols <> 1) || (y.Rows <> 1 && y.Cols <> 1) then raise (new Exception())
       let mutable ret : float = 0.0
       if x.Rows = 1 then
           if y.Rows = 1 then
               if x.Cols <> y.Cols then raise (new Exception())
               for i = 0 to x.Cols - 1 do ret <- ret + x.m.[0, i] * y.m.[0, i]
           else // y.Cols = 1
               if x.Cols <> y.Rows then raise (new Exception())
               for i = 0 to x.Cols - 1 do ret <- ret + x.m.[0, i] * y.m.[i, 0]
       else // x.Cols = 1
           if y.Rows = 1 then
               if x.Rows <> y.Cols then raise (new Exception())
               for i = 0 to x.Rows - 1 do ret <- ret + x.m.[i, 0] * y.m.[0, i]
           else // y.Cols = 1
               if x.Rows <> y.Rows then raise (new Exception())
               for i = 0 to x.Rows - 1 do ret <- ret + x.m.[i, 0] * y.m.[i, 0]
       ret
   member x.Norm () : float = sqrt(Matrix.Dot x x)
   member x.Normalize () : Matrix = x / x.Norm()
   member x.Each (f : float -> float) : Matrix =
       let ret = new Matrix(x.Rows, x.Cols)
       for i = 0 to x.Rows - 1 do
           for j = 0 to x.Cols - 1 do
               ret.m.[i, j] <- f x.m.[i, j]
       ret
end

let mat = new Matrix(2, 2)
mat.[1, 1] <- 3.0
let t = mat.[1, 1]
let mat' = mat.T()
let tmp = new Matrix(2, 3) * new Matrix(3, 4)
let tmp' = Matrix.Identity(3) * Matrix.Identity(3)
let tmp'' = Matrix.Dot (new Matrix(1, 3)) (new Matrix(3, 1))
let mat'' = mat.Each (fun x -> x + 1.0)