58 lines
1.1 KiB
Text
58 lines
1.1 KiB
Text
(* identity matrix *)
|
|
let eye n =
|
|
let a = Array.make_matrix n n 0.0 in
|
|
for i=0 to n-1 do
|
|
a.(i).(i) <- 1.0
|
|
done;
|
|
(a)
|
|
;;
|
|
|
|
(* matrix dimensions *)
|
|
let dim a = Array.length a, Array.length a.(0);;
|
|
|
|
(* make matrix from list in row-major order *)
|
|
let matrix p q v =
|
|
if (List.length v) <> (p * q)
|
|
then failwith "bad dimensions"
|
|
else
|
|
let a = Array.make_matrix p q (List.hd v) in
|
|
let rec g i j = function
|
|
| [] -> a
|
|
| x::v ->
|
|
a.(i).(j) <- x;
|
|
if j+1 < q
|
|
then g i (j+1) v
|
|
else g (i+1) 0 v
|
|
in
|
|
g 0 0 v
|
|
;;
|
|
|
|
(* matrix product *)
|
|
let matmul a b =
|
|
let n, p = dim a
|
|
and q, r = dim b in
|
|
if p <> q then failwith "bad dimensions" else
|
|
let c = Array.make_matrix n r 0.0 in
|
|
for i=0 to n-1 do
|
|
for j=0 to r-1 do
|
|
for k=0 to p-1 do
|
|
c.(i).(j) <- c.(i).(j) +. a.(i).(k) *. b.(k).(j)
|
|
done
|
|
done
|
|
done;
|
|
(c)
|
|
;;
|
|
|
|
(* generic exponentiation, usual algorithm *)
|
|
let pow one mul a n =
|
|
let rec g p x = function
|
|
| 0 -> x
|
|
| i ->
|
|
g (mul p p) (if i mod 2 = 1 then mul p x else x) (i/2)
|
|
in
|
|
g a one n
|
|
;;
|
|
|
|
(* example with integers *)
|
|
pow 1 ( * ) 2 16;;
|
|
(* - : int = 65536 *)
|