(* All the code with array input and sample generator *) module type DivAlgebra = sig type kind val sg : unit -> kind val conj : kind -> kind val zero : kind val one : kind val zeroQ : kind -> bool val (+) : kind -> kind -> kind val (-) : kind -> kind -> kind val ( * ) : kind -> kind -> kind val inv : kind -> kind val str : kind -> string val inp : float array -> kind end;; module RND = (struct (* Generate Normal Random Deviates *) let p = ref false and v = ref 0. and dr _ = (Random.float 2.) -. 1. let rec rnd d = if !p then (p := false; !v) else ( let x = dr () and y = dr () in let r2 = x *. x +. y *. y in if r2 < 1. then ( let q = sqrt (-2. *. (log r2) /. r2) in v := x *. q; p := true; y*.q) else rnd ()) end : sig val rnd : unit -> float end);; module BareReals = struct type kind = float let sg x = RND.rnd () let conj x = x let zero = 0. let one = 1. let zeroQ x = (abs_float x) < 1.e-10 let (+) = (+.) let (-) = (-.) let ( * ) = ( *.) let inv x = 1. /. x let str = string_of_float let inp x = x.(0) end;; (* And now, alternative functors *) module G = functor (Alg: DivAlgebra) -> struct open Alg type kind = {r: Alg.kind; i: Alg.kind} let sg x = {r = sg (); i = sg ()} and conj x = {r = conj x.r; i = zero - x.i} and zero = {r = zero; i = zero} and one = {r = one; i = zero} and zeroQ x = zeroQ x.r & zeroQ x.i and (+) x y = {r = x.r + y.r; i = x.i + y.i} and (-) x y = {r = x.r - y.r; i = x.i - y.i} and ( * ) x y = {r = x.r * y.r - (conj y.i) * x.i; i = y.i * x.r + x.i * (conj y.r)} and inv x = let d = inv (x.r * (conj x.r) + x.i * (conj x.i)) in {r = d * (conj x.r); i = zero - d * x.i} and str x = str x.r ^ ", " ^ str x.i and inp x = let hl = (Array.length x)/2 in {r = inp (Array.sub x 0 hl); i = inp (Array.sub x hl hl)} end;; module G = functor (Alg: DivAlgebra) -> struct open Alg type kind = (Alg.kind * Alg.kind) let sg x = sg (), sg () and conj (x, y) = conj x, zero - y and zero = zero, zero and one = one, zero and zeroQ (x, y) = zeroQ x & zeroQ y and (+) (a, b) (c, d) = a + c, b + d and (-) (a, b) (c, d) = a - c, b - d and ( * ) (a, b) (c, d) = a * c - (conj d) * b, d * a + b * (conj c) and inv (a, b) = let d = inv (a * (conj a) + b * (conj b)) in d * (conj a), zero - d * b and str (a, b) = str a ^ ", " ^ str b and inp x = let hl = (Array.length x)/2 in inp (Array.sub x 0 hl), inp (Array.sub x hl hl) end;; (* End of alternative functors *) module Reals = (BareReals : DivAlgebra);; module Quaternion = G(G(Reals));; Quaternion.str (Quaternion.inp [|0.; 1.; 2.; 3.|]);; (* yields - : string = "0., 1., 2., 3." *) let open Quaternion in let a = inp [|0.; 0.; 1.; 0.|] in str (a*a);; (* Yields - : string = "-1., 0., 0., 0." *) module Test = functor (Alg: DivAlgebra) -> struct open Alg let a = sg () and b = sg () and c = sg () let ai = inv a and ta value prop = (if (zeroQ value) then "" else "dont ") ^ prop let result = [ta ((a + (b + c)) - ((a + b) + c)) "+ associate"; ta ((conj a) * (conj b) - (conj (b * a))) "conjugate"; ta ((a * (b + c)) - (a * b + a * c)) "left distribute"; ta (((a + b) * c) - (a * c + b * c)) "right distribute"; ta ((b * a) * ai - b) "undo multiply"; ta (one - ai * a) "have * inverse"; ta (a * (b * c) - (a * b) * c) "* associate"; ta (a * b - b * a) "* commute"; ta ((a * a) * b - a * (a * b)) "left alternate"; ta ((a * b) * b - a * (b * b)) "right alternate"; ta (a * (ai * b) - b) "solve"] end;; module Tr = Test(Reals);; (* Reals *) Tr.result;; module Tr = Test(G(Reals));; (* Complexes *) Tr.result;; module Tr = Test(G(G(Reals)));; (* Quaternions *) Tr.result;; module Tr = Test(G(G(G(Reals))));; (* Octonions *) Tr.result;; module Tr = Test(G(G(G(G(Reals)))));; (* Sedenions *) Tr.result;; (* Two random normal deviate (RND) tests let sum = ref 0. and sq = ref 0. and aug x y = x := !x +. y in ( for j = 1 to 1000000 do let w = RND.rnd () in (aug sum w; aug sq (w*. w)) done; (!sum, !sq));; let h = Array.create 120 0 in ( for i = 1 to 1000000 do let x = 60 + int_of_float (floor (10. *. RND.rnd ())) in h.(x) <- h.(x) + 1 done; h);; *)