(* 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);;
*)