(* Obtained at http://www.arrakis.es/~worm/ *)

signature MONO_VECTOR =
  sig
    type vector
    type elem
    val maxLen : int
    val fromList : elem list -> vector
    val tabulate : (int * (int -> elem)) -> vector
    val length : vector -> int
    val sub : (vector * int) -> elem
    val extract : (vector * int * int option) -> vector
    val concat : vector list -> vector
    val mapi : ((int * elem) -> elem) -> (vector * int * int option) -> vector
    val map : (elem -> elem) -> vector -> vector
    val appi : ((int * elem) -> unit) -> (vector * int * int option) -> unit
    val app : (elem -> unit) -> vector -> unit
    val foldli : ((int * elem * 'a) -> 'a) -> 'a -> (vector * int * int option) -> 'a
    val foldri : ((int * elem * 'a) -> 'a) -> 'a -> (vector * int * int option) -> 'a
    val foldl : ((elem * 'a) -> 'a) -> 'a -> vector -> 'a
    val foldr : ((elem * 'a) -> 'a) -> 'a -> vector -> 'a 
  end

(*
 Copyright (c) Juan Jose Garcia Ripoll.
 All rights reserved.

 Refer to the COPYRIGHT file for license conditions
*)

(* COPYRIGHT
 
Redistribution and use in source and binary forms, with or
without modification, are permitted provided that the following
conditions are met:

1. Redistributions of source code must retain the above copyright
   notice, this list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above
   copyright notice, this list of conditions and the following
   disclaimer in the documentation and/or other materials provided
   with the distribution.

3. All advertising materials mentioning features or use of this
   software must display the following acknowledgement:
        This product includes software developed by Juan Jose
        Garcia Ripoll.

4. The name of Juan Jose Garcia Ripoll may not be used to endorse
   or promote products derived from this software without
   specific prior written permission.

THIS SOFTWARE IS PROVIDED BY JUAN JOSE GARCIA RIPOLL ``AS IS''
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL HE BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
OF SUCH DAMAGE.
*)

structure EvalTimer =
    struct
        local
            val TIME = ref (Time.now())
        in
            fun timerOn () =
                (TIME := Time.now(); ())
            fun timerRead () =
                Time.toMilliseconds(Time.-(Time.now(),!TIME))
            fun timerOff () =
                let val delta = timerRead()
                in
                    print "Elapsed: ";
                    print (LargeInt.toString delta);
                    print " ms\n"
                end
            fun time f = (timerOn(); f(); timerOff())
        end
    end
structure Loop =
    struct
        fun all (a, b, f) =
            if a > b then
                true
            else if f a then
                all (a+1, b, f)
            else
                false

        fun any (a, b, f) =
            if a > b then
                false
            else if f a then
                true
            else
                any (a+1, b, f)

        fun app (a, b, f) =
            if a < b then
                (f a; app (a+1, b, f))
            else
                ()

        fun app' (a, b, d, f) =
            if a < b then
                (f a; app' (a+d, b, d, f))
            else
                ()

        fun appi' (a, b, d, f) =
            if a < b then
                (f a; appi' (a+d, b, d, f))
            else
                ()
    end
(*
  INDEX         -Signature-

  Indices are a enumerable finite set of data with an order and a map
  to a continous nonnegative interval of integers.  In the sample
  implementation, Index, each index is a list of integers,
        [i1,...,in]
  and each set of indices is defined by a shape, which has the same
  shape of an index but with each integer incremented by one
        shape = [k1,...,kn]
        0 <= i1 < k1

  type storage = RowMajor | ColumnMajor
  order : storage
        Identifies:
                1) the underlying algorithms for this structure
                2) the most significant index
                3) the index that varies more slowly
                4) the total order
        RowMajor means that first index is most significant and varies
        more slowly, while ColumnMajor means that last index is the most
        significant and varies more slowly. For instance
                RowMajor => [0,0]<[0,1]<[1,0]<[1,1] (C, C++, Pascal)
                ColumnMajor => [0,0]>[1,0]>[0,1]>[1,1] (Fortran)
  last shape
  first shape
        Returns the last/first index that belongs to the sed defined by
        'shape'.
  inBounds shape index
        Checkes whether 'index' belongs to the set defined by 'shape'.
  toInt shape index
        As we said, indices can be sorted and mapped to a finite set of
        integers. 'toInt' obtaines the integer number that corresponds to
        a certain index.
  indexer shape
        It is equivalent to the partial evaluation 'toInt shape' but
        optimized for 'shape'.

  next shape index
  prev shape index
  next' shape index
  prev' shape index
        Obtain the following or previous index to the one we supply.
        next and prev return an object of type 'index option' so that
        if there is no such following/previous, the output is NONE.
        On the other hand, next'/prev' raise an exception when the
        output is not well defined and their output is always of type
        index. next/prev/next'/prev' raise an exception if 'index'
        does not belong to the set of 'shape'.

  all shape f
  any shape f
  app shape f
        Iterates 'f' over every index of the set defined by 'shape'.
        'all' stops when 'f' first returns false, 'any' stops when
        'f' first returns true and 'app' does not stop and discards the
        output of 'f'.

  compare(a,b)
        Returns LESS/GREATER/EQUAL according to the total order which
        is defined in the set of all indices.
  <,>,eq,<=,>=,<>
        Reduced comparisons which are defined in terms of 'compare'.

  validShape t
  validIndex t
        Checks whether 't' conforms a valid shape or index.

  iteri shape f
*)

signature INDEX =
    sig
        type t
        type indexer = t -> int
        datatype storage = RowMajor | ColumnMajor

        exception Index
        exception Shape

        val order : storage
        val toInt : t -> t -> int
        val length : t -> int
        val first : t -> t
        val last : t -> t
        val next : t -> t -> t option
        val prev : t -> t -> t option
        val next' : t -> t -> t
        val prev' : t -> t -> t
        val indexer : t -> (t -> int)

        val inBounds : t -> t -> bool
        val compare : t * t -> order
        val < : t * t -> bool
        val > : t * t -> bool
        val eq : t * t -> bool
        val <= : t * t -> bool
        val >= : t * t -> bool
        val <> : t * t -> bool
        val - : t * t -> t

        val validShape : t -> bool
        val validIndex : t -> bool

        val all : t -> (t -> bool) -> bool
        val any : t -> (t -> bool) -> bool
        val app : t -> (t -> unit) -> unit
    end
structure Index : INDEX =
    struct
        type t = int list
        type indexer = t -> int
        datatype storage = RowMajor | ColumnMajor

        exception Index
        exception Shape

        val order = ColumnMajor

        fun validShape shape = List.all (fn x => x > 0) shape

        fun validIndex index = List.all (fn x => x >= 0) index

        fun toInt shape index =
            let fun loop ([], [], accum, _) = accum
                  | loop ([], _, _, _) = raise Index
                  | loop (_, [], _, _) = raise Index
                  | loop (i::ri, l::rl, accum, fac) =
                if (i >= 0) andalso (i < l) then
                    loop (ri, rl, i*fac + accum, fac*l)
                else
                    raise Index
            in loop (index, shape, 0, 1)
            end

        (* ----- CACHED LINEAR INDEXER -----

           An indexer is a function that takes a list of
           indices, validates it and produces a nonnegative
           integer number. In short, the indexer is the
           mapper from indices to element positions in
           arrays.

           'indexer' builds such a mapper by optimizing
           the most common cases, which are 1d and 2d
           tensors.
         *)
    local
        fun doindexer [] _ = raise Shape
          | doindexer [a] [dx] =
            let fun f [x] = if (x > 0) andalso (x < a)
                            then x
                            else raise Index
                  | f _ = raise Index
            in f end
          | doindexer [a,b] [dx, dy] =
            let fun f [x,y] = if ((x > 0) andalso (x < a) andalso
                                  (y > 0) andalso (y < b))
                              then x + dy * y
                              else raise Index
                  | f _ = raise Index
            in f end
          | doindexer [a,b,c] [dx,dy,dz] =
            let fun f [x,y,z] = if ((x > 0) andalso (x < a) andalso
                                    (y > 0) andalso (y < b) andalso
                                    (z > 0) andalso (z < c))
                                then x + dy * y + dz * z
                                else raise Index
                  | f _ = raise Index
            in f end
          | doindexer shape memo =
            let fun f [] [] accum [] = accum
                  | f _  _  _ [] = raise Index
                  | f (fac::rf) (ndx::ri) accum (dim::rd) =
                    if (ndx >= 0) andalso (ndx < dim) then
                        f rf ri (accum + ndx * fac) rd
                    else
                        raise Index
            in f shape memo 0
            end
    in
        fun indexer shape =
            let fun memoize accum [] = []
                  | memoize accum (dim::rd) =
                    accum :: (memoize (dim * accum) rd)
            in
                if validShape shape
                then doindexer shape (memoize 1 shape)
                else raise Shape
            end
    end

        fun length shape =
            let fun prod (a,b) =
                if b < 0 then raise Shape else a * b
            in foldl prod 1 shape
            end

        fun first shape = map (fn x => 0) shape

        fun last [] = []
          | last (size :: rest) =
            if size < 1
            then raise Shape
            else size - 1 :: last rest

        fun next' [] [] = raise Subscript
          | next' _ [] = raise Index
          | next' [] _ = raise Index
          | next' (dimension::restd) (index::resti) =
            if (index + 1) < dimension
            then (index + 1) :: resti
            else 0 :: (next' restd resti)

        fun prev' [] [] = raise Subscript
          | prev' _ [] = raise Index
          | prev' [] _ = raise Index
          | prev' (dimension::restd) (index::resti) =
            if (index > 0)
            then index - 1 :: resti
            else dimension - 1 :: prev' restd resti

        fun next shape index = (SOME (next' shape index)) handle
            Subscript => NONE

        fun prev shape index = (SOME (prev' shape index)) handle
            Subscript => NONE

        fun inBounds shape index =
            ListPair.all (fn (x,y) => (x >= 0) andalso (x < y))
            (index, shape)

        fun compare ([],[]) = EQUAL
          | compare (_, []) = raise Index
          | compare ([],_) = raise Index
          | compare (a::ra, b::rb) =
            case Int.compare (a,b) of
                EQUAL => compare (ra,rb)
              | LESS => LESS
              | GREATER => GREATER

    local
        fun iterator a inner =
            let fun loop accum f =
                let fun innerloop i =
                    if i < a
                    then if inner (i::accum) f
                         then innerloop (i+1)
                         else false
                    else true
                in innerloop 0
                end
            in loop
            end
        fun build_iterator [a] =
            let fun loop accum f =
                let fun innerloop i =
                    if i < a
                    then if f (i::accum)
                         then innerloop (i+1)
                         else false
                    else true
                in innerloop 0
                end
            in loop
            end
          | build_iterator (a::rest) = iterator a (build_iterator rest)
    in
        fun all shape = build_iterator shape []
    end

    local
        fun iterator a inner =
            let fun loop accum f =
                let fun innerloop i =
                    if i < a
                    then if inner (i::accum) f
                         then true
                         else innerloop (i+1)
                    else false
                in innerloop 0
                end
            in loop
            end
        fun build_iterator [a] =
            let fun loop accum f =
                let fun innerloop i =
                    if i < a
                    then if f (i::accum)
                         then true
                         else innerloop (i+1)
                    else false
                in innerloop 0
                end
            in loop
            end
          | build_iterator (a::rest) = iterator a (build_iterator rest)
    in
        fun any shape = build_iterator shape []
    end

    local
        fun iterator a inner =
            let fun loop accum f =
                let fun innerloop i =
                    if i < a
                    then (inner (i::accum) f;
                          innerloop (i+1))
                    else ()
                in innerloop 0
                end
            in loop
            end
        fun build_iterator [a] =
            let fun loop accum f =
                let fun innerloop i =
                    if i < a
                    then (f (i::accum); innerloop (i+1))
                    else ()
                in innerloop 0
                end
            in loop
            end
          | build_iterator (a::rest) = iterator a (build_iterator rest)
    in
        fun app shape = build_iterator shape []
    end

        fun a < b = compare(a,b) = LESS
        fun a > b = compare(a,b) = GREATER
        fun eq (a, b) = compare(a,b) = EQUAL
        fun a <> b = not (a = b)
        fun a <= b = not (a > b)
        fun a >= b = not (a < b)
        fun a - b = ListPair.map Int.- (a,b)

    end
(*
 Copyright (c) Juan Jose Garcia Ripoll.
 All rights reserved.
 
 Refer to the COPYRIGHT file for license conditions
*)

(*
 TENSOR         - Signature -

 Polymorphic tensors of any type. With 'tensor' we denote a (mutable)
 array of any rank, with as many indices as one wishes, and that may
 be traversed (map, fold, etc) according to any of those indices.

 type 'a tensor
        Polymorphic tensor whose elements are all of type 'a.
 val storage = RowMajor | ColumnMajor
        RowMajor = data is stored in consecutive cells, first index
        varying fastest (FORTRAN convention)
        ColumnMajor = data is stored in consecutive cells, last
        index varying fastest (C,C++,Pascal,CommonLisp convention)
 new ([i1,...,in],init)
        Build a new tensor with n indices, each of sizes i1...in,
        filled with 'init'.
 fromArray (shape,data)
 fromList (shape,data)
        Use 'data' to fill a tensor of that shape. An exception is
        raised if 'data' is too large or too small to properly
        fill the vector. Later use of a 'data' array is disregarded
        -- one must think that the tensor now owns the array.
 length tensor
 rank tensor
 shape tensor
        Return the number of elements, the number of indices and
        the shape (size of each index) of the tensor.
 toArray tensor
        Return the data of the tensor in the form of an array.
        Mutation of this array may lead to unexpected behavior.

 sub (tensor,[i1,...,in])
 update (tensor,[i1,...,in],new_value)
        Access the element that is indexed by the numbers [i1,..,in]

 app f a
 appi f a
        The same as 'map' and 'mapi' but the function 'f' outputs
        nothing and no new array is produced, i.e. one only seeks
        the side effect that 'f' may produce.
 map2 operation a b
        Apply function 'f' to pairs of elements of 'a' and 'b'
        and build a new tensor with the output. Both operands
        must have the same shape or an exception is raised.
        The procedure is sequential, as specified by 'storage'.
 foldl operation a n
        Fold-left the elements of tensor 'a' along the n-th
        index.
 all test a
 any test a
        Folded boolean tests on the elements of the tensor.
*)

signature TENSOR =
    sig
        structure Array : ARRAY
        structure Index : INDEX
        type index = Index.t
        type 'a tensor

        val new : index * 'a -> 'a tensor
        val tabulate : index * (index -> 'a) -> 'a tensor
        val length : 'a tensor -> int
        val rank : 'a tensor -> int
        val shape : 'a tensor -> (index)
        val reshape : index -> 'a tensor -> 'a tensor
        val fromList : index * 'a list -> 'a tensor
        val fromArray : index * 'a array -> 'a tensor
        val toArray : 'a tensor -> 'a array

        val sub : 'a tensor * index -> 'a
        val update : 'a tensor * index * 'a -> unit
        val map : ('a -> 'b) -> 'a tensor -> 'b tensor
        val map2 : ('a * 'b -> 'c) -> 'a tensor -> 'b tensor -> 'c tensor
        val app : ('a -> unit) -> 'a tensor -> unit
        val appi : (int * 'a -> unit) -> 'a tensor -> unit
        val foldl : ('c * 'a -> 'c) -> 'c -> 'a tensor -> int -> 'c tensor
        val all : ('a -> bool) -> 'a tensor -> bool
        val any : ('a -> bool) -> 'a tensor -> bool
    end

(*
 Copyright (c) Juan Jose Garcia Ripoll.
 All rights reserved.
 
 Refer to the COPYRIGHT file for license conditions
*)

structure Tensor : TENSOR =
    struct
        structure Array = Array
        structure Index = Index
            
        type index = Index.t
        type 'a tensor = {shape : index, indexer : Index.indexer, data : 'a array}

        exception Shape
        exception Match
        exception Index

    local
    (*----- LOCALS -----*)

        fun make' (shape, data) =
            {shape = shape, indexer = Index.indexer shape, data = data}

        fun toInt {shape, indexer, data} index = indexer index

        fun array_map f a =
            let fun apply index = f(Array.sub(a,index)) in
                Array.tabulate(Array.length a, apply)
            end

        fun splitList (l as (a::rest), place) =
            let fun loop (left,here,right) 0 =  (List.rev left,here,right)
                  | loop (_,_,[]) place = raise Index
                  | loop (left,here,a::right) place = 
                loop (here::left,a,right) (place-1)
            in
                if place <= 0 then
                    loop ([],a,rest) (List.length rest - place)
                else
                    loop ([],a,rest) (place - 1)
            end

    in
    (*----- STRUCTURAL OPERATIONS & QUERIES ------*)

        fun new (shape, init) =
            if not (Index.validShape shape) then
                raise Shape
            else
                let val length = Index.length shape in
                    {shape = shape,
                     indexer = Index.indexer shape,
                     data = Array.array(length,init)}
                end

        fun toArray {shape, indexer, data} = data

        fun length {shape, indexer, data} =  Array.length data

        fun shape {shape, indexer, data} = shape

        fun rank t = List.length (shape t)

        fun reshape new_shape tensor =
            if Index.validShape new_shape then
                case (Index.length new_shape) = length tensor of
                    true => make'(new_shape, toArray tensor)
                  | false => raise Match
            else
                raise Shape

        fun fromArray (s, a) =
            case Index.validShape s andalso 
                 ((Index.length s) = (Array.length a)) of
                 true => make'(s, a)
               | false => raise Shape

        fun fromList (s, a) = fromArray (s, Array.fromList a)

        fun tabulate (shape,f) =
            if Index.validShape shape then
                let val last = Index.last shape
                    val length = Index.length shape
                    val c = Array.array(length, f last)
                    fun dotable (c, indices, i) =
                        (Array.update(c, i, f indices);
                         case i of
                             0 => c
                           | i => dotable(c, Index.prev' shape indices, i-1))
                in
                    make'(shape,dotable(c, Index.prev' shape last, length-1))
                end
            else
                raise Shape

        (*----- ELEMENTWISE OPERATIONS -----*)

        fun sub (t, index) = Array.sub(#data t, toInt t index)

        fun update (t, index, value) =
            Array.update(toArray t, toInt t index, value)

        fun map f {shape, indexer, data} =
            {shape = shape, indexer = indexer, data = array_map f data}

        fun map2 f t1 t2=
            let val {shape, indexer, data} = t1
                val {shape=shape2, indexer=indexer2, data=data2} = t2
                fun apply i = f (Array.sub(data,i), Array.sub(data2,i))
                val len = Array.length data
            in
                if Index.eq(shape, shape2) then
                    {shape = shape,
                     indexer = indexer,
                     data = Array.tabulate(len, apply)}
                else
                    raise Match
        end

        fun appi f tensor = Array.appi f (toArray tensor)

        fun app f tensor = Array.app f (toArray tensor)

        fun all f tensor =
            let val a = toArray tensor
            in Loop.all(0, length tensor - 1, fn i =>
                        f (Array.sub(a, i)))
            end

        fun any f tensor =
            let val a = toArray tensor
            in Loop.any(0, length tensor - 1, fn i =>
                        f (Array.sub(a, i)))
            end

        fun foldl f init {shape, indexer, data=a} index =
            let val (head,lk,tail) = splitList(shape, index)
                val li = Index.length head
                val lj = Index.length tail
                val c = Array.array(li * lj,init)
                fun loopi (0, _,  _)  = ()
                  | loopi (i, ia, ic) =
                    (Array.update(c, ic, f(Array.sub(c,ic), Array.sub(a,ia)));
                     loopi (i-1, ia+1, ic+1))
                fun loopk (0, ia, _)  = ia
                  | loopk (k, ia, ic) = (loopi (li, ia, ic);
                                         loopk (k-1, ia+li, ic))
                fun loopj (0, _,  _)  = ()
                  | loopj (j, ia, ic) = loopj (j-1, loopk(lk,ia,ic), ic+li)
            in
                loopj (lj, 0, 0);
                make'(head @ tail, c)
            end

    end
    end (* Tensor *)

(*
 Copyright (c) Juan Jose Garcia Ripoll.
 All rights reserved.
 
 Refer to the COPYRIGHT file for license conditions
*)

(*
 MONO_TENSOR            - signature -

 Monomorphic tensor of arbitrary data (not only numbers). Operations
 should be provided to run the data in several ways, according to one
 index.

 type tensor
        The type of the tensor itself
 type elem
        The type of every element
 val storage = RowMajor | ColumnMajor
        RowMajor = data is stored in consecutive cells, first index
        varying fastest (FORTRAN convention)
        ColumnMajor = data is stored in consecutive cells, last
        index varying fastest (C,C++,Pascal,CommonLisp convention)
 new ([i1,...,in],init)
        Build a new tensor with n indices, each of sizes i1...in,
        filled with 'init'.
 fromArray (shape,data)
 fromList (shape,data)
        Use 'data' to fill a tensor of that shape. An exception is
        raised if 'data' is too large or too small to properly
        fill the vector. Later use of a 'data' array is disregarded
        -- one must think that the tensor now owns the array.
 length tensor
 rank tensor
 shape tensor
        Return the number of elements, the number of indices and
        the shape (size of each index) of the tensor.
 toArray tensor
        Return the data of the tensor in the form of an array.
        Mutation of this array may lead to unexpected behavior.
        The data in the array is stored according to `storage'.

 sub (tensor,[i1,...,in])
 update (tensor,[i1,...,in],new_value)
        Access the element that is indexed by the numbers [i1,..,in]

 map f a
 mapi f a
        Produce a new array by mapping the function sequentially
        as specified by 'storage', to each element of tensor 'a'.
        In 'mapi' the function receives a (indices,value) tuple,
        while in 'map' it only receives the value.
 app f a
 appi f a
        The same as 'map' and 'mapi' but the function 'f' outputs
        nothing and no new array is produced, i.e. one only seeks
        the side effect that 'f' may produce.
 map2 operation a b
        Apply function 'f' to pairs of elements of 'a' and 'b'
        and build a new tensor with the output. Both operands
        must have the same shape or an exception is raised.
        The procedure is sequential, as specified by 'storage'.
 foldl operation a n
        Fold-left the elements of tensor 'a' along the n-th
        index.
 all test a
 any test a
        Folded boolean tests on the elements of the tensor.

 map', map2', foldl'
        Polymorphic versions of map, map2, foldl.
*)

signature MONO_TENSOR =
    sig
        structure Array : MONO_ARRAY
        structure Index : INDEX
        type index = Index.t
        type elem
        type tensor
        type t = tensor

        val new : index * elem -> tensor
        val tabulate : index * (index -> elem) -> tensor
        val length : tensor -> int
        val rank : tensor -> int
        val shape : tensor -> (index)
        val reshape : index -> tensor -> tensor
        val fromList : index * elem list -> tensor
        val fromArray : index * Array.array -> tensor
        val toArray : tensor -> Array.array

        val sub : tensor * index -> elem
        val update : tensor * index * elem -> unit
        val map : (elem -> elem) -> tensor -> tensor
        val map2 : (elem * elem -> elem) -> tensor -> tensor -> tensor
        val app : (elem -> unit) -> tensor -> unit
        val appi : (int * elem -> unit) -> tensor -> unit
        val foldl : (elem * 'a -> 'a) -> 'a -> tensor -> tensor
        val foldln : (elem * elem -> elem) -> elem -> tensor -> int -> tensor
        val all : (elem -> bool) -> tensor -> bool
        val any : (elem -> bool) -> tensor -> bool

        val map' : (elem -> 'a) -> tensor -> 'a Tensor.tensor
        val map2' : (elem * elem -> 'a) -> tensor -> tensor -> 'a Tensor.tensor
        val foldl' : ('a * elem -> 'a) -> 'a -> tensor -> int -> 'a Tensor.tensor
    end

(*
 NUMBER         - Signature -

 Guarantees a structure with a minimal number of mathematical operations
 so as to build an algebraic structure named Tensor.
 *)

signature NUMBER =
    sig
        type t
        val zero : t
        val one  : t
        val ~ : t -> t
        val + : t * t -> t
        val - : t * t -> t
        val * : t * t -> t
        val / : t * t -> t
        val toString : t -> string
    end

signature NUMBER =
    sig
        type t
        val zero : t
        val one : t

        val + : t * t -> t
        val - : t * t -> t
        val * : t * t -> t
        val *+ : t * t * t -> t
        val *- : t * t * t -> t
        val ** : t * int -> t

        val ~ : t -> t
        val abs : t -> t
        val signum : t -> t

        val == : t * t -> bool
        val != : t * t -> bool

        val toString : t -> string
        val fromInt : int -> t
        val scan : (char,'a) StringCvt.reader -> (t,'a) StringCvt.reader
    end

signature INTEGRAL_NUMBER =
    sig
        include NUMBER

        val quot : t * t -> t
        val rem  : t * t -> t
        val mod  : t * t -> t
        val div  : t * t -> t

        val compare : t * t -> order
        val < : t * t -> bool
        val > : t * t -> bool
        val <= : t * t -> bool
        val >= : t * t -> bool

        val max : t * t -> t
        val min : t * t -> t
    end

signature FRACTIONAL_NUMBER =
    sig
        include NUMBER

        val pi : t
        val e : t

        val / : t * t -> t
        val recip : t -> t

        val ln : t -> t
        val pow : t * t -> t
        val exp : t -> t
        val sqrt : t -> t

        val cos : t -> t
        val sin : t -> t
        val tan : t -> t
        val sinh : t -> t
        val cosh : t -> t
        val tanh : t -> t

        val acos : t -> t
        val asin : t -> t
        val atan : t -> t
        val asinh : t -> t
        val acosh : t -> t
        val atanh : t -> t
        val atan2 : t * t -> t
    end

signature REAL_NUMBER =
    sig
        include FRACTIONAL_NUMBER

        val compare : t * t -> order
        val < : t * t -> bool
        val > : t * t -> bool
        val <= : t * t -> bool
        val >= : t * t -> bool

        val max : t * t -> t
        val min : t * t -> t
    end

signature COMPLEX_NUMBER =
    sig
        include FRACTIONAL_NUMBER

        structure Real : REAL_NUMBER
        type real = Real.t

        val make : real * real -> t
        val split : t -> real * real
        val realPart : t -> real
        val imagPart : t -> real
        val abs2 : t -> real
    end

structure INumber : INTEGRAL_NUMBER =
    struct
        open Int
        type t = Int.int
        val zero = 0
        val one = 1

        infix **
        fun i ** n =
            let fun loop 0 = 1
                  | loop 1 = i
                  | loop n =
                let val x = loop (Int.div(n, 2))
                    val m = Int.mod(n, 2)
                in
                    if m = 0 then
                        x * x
                    else
                        x * x * i
                end
            in if n < 0
               then raise Domain
               else loop n
            end

        fun signum i = case compare(i, 0) of
            GREATER => 1
          | EQUAL => 0
          | LESS => ~1

        infix ==
        infix !=
        fun a == b = a = b
        fun a != b = (a <> b)
        fun *+(b,c,a) = b * c + a
        fun *-(b,c,a) = b * c - b

        fun scan getc = Int.scan StringCvt.DEC getc
    end

structure RNumber : REAL_NUMBER =
    struct
        open Real
        open Real.Math
        type t = Real.real
        val zero = 0.0
        val one = 1.0

        fun signum x = case compare(x,0.0) of
            LESS => ~1.0
          | GREATER => 1.0
          | EQUAL => 0.0

        fun recip x = 1.0 / x

        infix **
        fun i ** n =
            let fun loop 0 = one
                  | loop 1 = i
                  | loop n =
                let val x = loop (Int.div(n, 2))
                    val m = Int.mod(n, 2)
                in
                    if m = 0 then
                        x * x
                    else
                        x * x * i
                end
            in if Int.<(n, 0)
               then raise Domain
               else loop n
            end

        fun max (a, b) = if a < b then b else a
        fun min (a, b) = if a < b then a else b

        fun asinh x = ln (x + sqrt(1.0 + x * x))
        fun acosh x = ln (x + (x + 1.0) * sqrt((x - 1.0)/(x + 1.0)))
        fun atanh x = ln ((1.0 + x) / sqrt(1.0 - x * x))

    end
(*
 Complex(R)     - Functor -

 Provides support for complex numbers based on tuples. Should be
 highly efficient as most operations can be inlined.
 *)

structure CNumber : COMPLEX_NUMBER =
struct
        structure Real = RNumber

        type t = Real.t * Real.t
        type real = Real.t

        val zero = (0.0,0.0)
        val one = (1.0,0.0)
        val pi = (Real.pi, 0.0)
        val e = (Real.e, 0.0)

        fun make (r,i) = (r,i) : t
        fun split z = z
        fun realPart (r,_) = r
        fun imagPart (_,i) = i

        fun abs2 (r,i) = Real.+(Real.*(r,r),Real.*(i,i)) (* FIXME!!! *)
        fun arg (r,i) = Real.atan2(i,r)
        fun modulus z = Real.sqrt(abs2 z)
        fun abs z = (modulus z, 0.0)
        fun signum (z as (r,i)) =
            let val m = modulus z
            in (Real./(r,m), Real./(i,m))
            end

        fun ~ (r1,i1) = (Real.~ r1, Real.~ i1)
        fun (r1,i1) + (r2,i2) = (Real.+(r1,r2), Real.+(i1,i2))
        fun (r1,i1) - (r2,i2) = (Real.-(r1,r2), Real.-(i1,i1))
        fun (r1,i1) * (r2,i2) = (Real.-(Real.*(r1,r2),Real.*(i1,i2)),
                                 Real.+(Real.*(r1,i2),Real.*(r2,i1)))
        fun (r1,i1) / (r2,i2) =
            let val modulus = abs2(r2,i2)
                val (nr,ni) = (r1,i1) * (r2,i2)
            in
                (Real./(nr,modulus), Real./(ni,modulus))
            end
        fun *+((r1,i1),(r2,i2),(r0,i0)) =
            (Real.*+(Real.~ i1, i2, Real.*+(r1,r2,r0)),
             Real.*+(r2, i2, Real.*+(r1,i2,i0)))
        fun *-((r1,i1),(r2,i2),(r0,i0)) =
            (Real.*+(Real.~ i1, i2, Real.*-(r1,r2,r0)),
             Real.*+(r2, i2, Real.*-(r1,i2,i0)))

        infix **
        fun i ** n =
            let fun loop 0 = one
                  | loop 1 = i
                  | loop n =
                let val x = loop (Int.div(n, 2))
                    val m = Int.mod(n, 2)
                in
                    if m = 0 then
                        x * x
                    else
                        x * x * i
                end
            in if Int.<(n, 0)
                   then raise Domain
               else loop n
            end

        fun recip (r1, i1) = 
            let val modulus = abs2(r1, i1)
            in (Real./(r1, modulus), Real./(Real.~ i1, modulus))
            end
        fun ==(z, w) = Real.==(realPart z, realPart w) andalso Real.==(imagPart z, imagPart w)
        fun !=(z, w) = Real.!=(realPart z, realPart w) andalso Real.!=(imagPart z, imagPart w)
        fun fromInt i = (Real.fromInt i, 0.0)
        fun toString (r,i) =
            String.concat ["(",Real.toString r,",",Real.toString i,")"]

        fun exp (x, y) =
            let val expx = Real.exp x
            in (Real.*(x, (Real.cos y)), Real.*(x, (Real.sin y)))
            end

    local
        val half = Real.recip (Real.fromInt 2)
    in
        fun sqrt (z as (x,y)) =
            if Real.==(x, 0.0) andalso Real.==(y, 0.0) then
                zero
            else
                let val m = Real.+(modulus z, Real.abs x)
                    val u' = Real.sqrt (Real.*(m, half))
                    val v' = Real./(Real.abs y , Real.+(u',u'))
                    val (u,v) = if Real.<(x, 0.0) then (v',u') else (u',v')
                in (u, if Real.<(y, 0.0) then Real.~ v else v)
                end
    end
        fun ln z = (Real.ln (modulus z), arg z)

        fun pow (z, n) =
            let val l = ln z
            in exp (l * n)
            end

        fun sin (x, y) = (Real.*(Real.sin x, Real.cosh y),
                          Real.*(Real.cos x, Real.sinh y))
        fun cos (x, y) = (Real.*(Real.cos x, Real.cosh y),
                          Real.~ (Real.*(Real.sin x, Real.sinh y)))
        fun tan (x, y) =
            let val (sx, cx) = (Real.sin x, Real.cos x)
                val (shy, chy) = (Real.sinh y, Real.cosh y)
                val a = (Real.*(sx, chy), Real.*(cx, shy))
                val b = (Real.*(cx, chy), Real.*(Real.~ sx, shy))
            in a / b
            end

        fun sinh (x, y) = (Real.*(Real.cos y, Real.sinh x),
                           Real.*(Real.sin y, Real.cosh x))
        fun cosh (x, y) = (Real.*(Real.cos y, Real.cosh x),
                           Real.*(Real.sin y, Real.sinh x))
        fun tanh (x, y) =
            let val (sy, cy) = (Real.sin y, Real.cos y)
                val (shx, chx) = (Real.sinh x, Real.cosh x)
                val a = (Real.*(cy, shx), Real.*(sy, chx))
                val b = (Real.*(cy, chx), Real.*(sy, shx))
            in a / b
            end

        fun asin (z as (x,y)) =
            let val w = sqrt (one - z * z)
                val (x',y') = ln ((Real.~ y, x) + w)
            in (y', Real.~ x')
            end

        fun acos (z as (x,y)) = 
            let val (x', y') = sqrt (one + z * z)
                val (x'', y'') = ln (z + (Real.~ y', x'))
            in (y'', Real.~ x'')
            end

        fun atan (z as (x,y)) =
            let val w = sqrt (one + z*z)
                val (x',y') = ln ((Real.-(1.0, y), x) / w)
            in (y', Real.~ x')
            end

        fun atan2 (y, x) = atan(y / x)

        fun asinh x = ln (x + sqrt(one + x * x))
        fun acosh x = ln (x + (x + one) * sqrt((x - one)/(x + one)))
        fun atanh x = ln ((one + x) / sqrt(one - x * x))

        fun scan getc =
            let val scanner = Real.scan getc
            in fn stream => 
                  case scanner stream of
                      NONE => NONE
                    | SOME (a, rest) =>
                      case scanner rest of
                          NONE => NONE
                        | SOME (b, rest) => SOME (make(a,b), rest)
            end

end (* ComplexNumber *)

(*
 Copyright (c) Juan Jose Garcia Ripoll.
 All rights reserved.
 
 Refer to the COPYRIGHT file for license conditions
*)

structure INumberArray =
    struct
        open Array
        type array = INumber.t array
        type vector = INumber.t vector
        type elem  = INumber.t
        structure Vector =
            struct
                open Vector
                type vector = INumber.t Vector.vector
                type elem = INumber.t
            end
        fun map f a = tabulate(length a, fn x => (f (sub(a,x))))
        fun mapi f a = tabulate(length a, fn x => (f (x,sub(a,x))))
        fun map2 f a b = tabulate(length a, fn x => (f(sub(a,x),sub(b,x))))
    end

structure RNumberArray =
    struct
        open Real64Array
        val sub = Unsafe.Real64Array.sub
        val update = Unsafe.Real64Array.update
        fun map f a = tabulate(length a, fn x => (f (sub(a,x))))
        fun mapi f a = tabulate(length a, fn x => (f (x,sub(a,x))))
        fun map2 f a b = tabulate(length a, fn x => (f(sub(a,x),sub(b,x))))
    end

(*--------------------- COMPLEX ARRAY -------------------------*)

structure BasicCNumberArray =
struct
        structure Complex : COMPLEX_NUMBER = CNumber
        structure Array : MONO_ARRAY = RNumberArray

        type elem = Complex.t
        type array = Array.array * Array.array

        val maxLen = Array.maxLen

        fun length (a,b) = Array.length a

        fun sub ((a,b),index) = Complex.make(Array.sub(a,index),Array.sub(b,index))

        fun update ((a,b),index,z) =
            let val (re,im) = Complex.split z in
                Array.update(a, index, re);
                Array.update(b, index, im)
            end

    local
        fun makeRange (a, start, NONE) = makeRange(a, start, SOME (length a - 1))
          | makeRange (a, start, SOME last) =
            let val len = length a
                val diff = last - start
            in
                if (start >= len) orelse (last >= len) then
                    raise Subscript
                else if diff < 0 then
                    (a, start, 0)
                else
                    (a, start, diff + 1)
            end

    in

        fun array (size,z:elem) =
            let val realsize = size * 2
                val r = Complex.realPart z
                val i = Complex.imagPart z in
                    (Array.array(size,r), Array.array(size,i))
            end

        fun zeroarray size =
            (Array.array(size,Complex.Real.zero),
             Array.array(size,Complex.Real.zero))

        fun tabulate (size,f) =
            let val a = array(size, Complex.zero)
                fun loop i =
                    case i = size of
                        true => a
                      | false => (update(a, i, f i); loop (i+1))
            in
                loop 0
            end

        fun fromList list =
            let val length = List.length list
                val a = zeroarray length
                fun loop (_, []) = a
                  | loop (i, z::rest) = (update(a, i, z);
                                         loop (i+1, rest))
            in
                loop(0,list)
            end

        fun extract range =
            let val (a, start, len) = makeRange range
                fun copy i = sub(a, i + start)
            in tabulate(len, copy)
            end

        fun concat array_list =
            let val total_length = foldl (op +) 0 (map length array_list)
                val a = array(total_length, Complex.zero)
                fun copy (_, []) = a
                  | copy (pos, v::rest) =
                    let fun loop i =
                        case i = 0 of
                            true => ()
                          | false => (update(a, i+pos, sub(v, i)); loop (i-1))
                    in (loop (length v - 1); copy(length v + pos, rest))
                    end
            in
                copy(0, array_list)
            end

        fun copy {src : array, si : int, len : int option, dst : array, di : int } =
            let val (a, ia, la) = makeRange (src, si, len)
                val (b, ib, lb) = makeRange (dst, di, len)
                fun copy i =
                    case i < 0 of
                        true => ()
                      | false => (update(b, i+ib, sub(a, i+ia)); copy (i-1))
            in copy (la - 1)
            end

        val copyVec = copy

        fun modifyi f range =
            let val (a, start, len) = makeRange range
                val last = start + len
                fun loop i =
                    case i >= last of
                        true => ()
                      | false => (update(a, i, f(i, sub(a,i))); loop (i+1))
            in loop start
            end

        fun modify f a =
            let val last = length a
                fun loop i =
                    case i >= last of
                        true => ()
                      | false => (update(a, i, f(sub(a,i))); loop (i+1))
            in loop 0
            end

        fun app f a =
            let val size = length a
                fun loop i =
                    case i = size of
                        true => ()
                      | false => (f(sub(a,i)); loop (i+1))
            in
                loop 0
            end

        fun appi f range =
            let val (a, start, len) = makeRange range
                val last = start + len
                fun loop i =
                    case i >= last of
                        true => ()
                      | false => (f(i, sub(a,i)); loop (i+1))
            in
                loop start
            end

        fun map f a =
            let val len = length a
                val c = zeroarray len
                fun loop ~1 = c
                  | loop i = (update(a, i, f(sub(a,i))); loop (i-1))
            in loop (len-1)
            end

        fun map2 f a b =
            let val len = length a
                val c = zeroarray len
                fun loop ~1 = c
                  | loop i = (update(c, i, f(sub(a,i),sub(b,i)));
                              loop (i-1))
            in loop (len-1)
            end

        fun mapi f range =
            let val (a, start, len) = makeRange range
                fun rule i = f (i+start, sub(a, i+start))
            in tabulate(len, rule)
            end

        fun foldli f init range =
            let val (a, start, len) = makeRange range
                val last = start + len - 1
                fun loop (i, accum) =
                    case i > last of
                        true => accum
                      | false => loop (i+1, f(i, sub(a,i), accum))
            in loop (start, init)
            end

        fun foldri f init range =
            let val (a, start, len) = makeRange range
                val last = start + len - 1
                fun loop (i, accum) =
                    case i < start of
                        true => accum
                      | false => loop (i-1, f(i, sub(a,i), accum))
            in loop (last, init)
            end

        fun foldl f init a = foldli (fn (_, a, x) => f(a,x)) init (a,0,NONE)
        fun foldr f init a = foldri (fn (_, x, a) => f(x,a)) init (a,0,NONE)
    end
end (* BasicCNumberArray *)


structure CNumberArray =
    struct
        structure Vector =
            struct
                open BasicCNumberArray
                type vector = array
            end : MONO_VECTOR
        type vector = Vector.vector
        open BasicCNumberArray
    end (* CNumberArray *)
structure INumber : INTEGRAL_NUMBER =
    struct
        open Int
        type t = Int.int
        val zero = 0
        val one = 1
        infix **
        fun i ** n =
            let fun loop 0 = 1
                  | loop 1 = i
                  | loop n =
                let val x = loop (Int.div(n, 2))
                    val m = Int.mod(n, 2)
                in
                    if m = 0 then
                        x * x
                    else
                        x * x * i
                end
            in if n < 0
               then raise Domain
               else loop n
            end
        fun signum i = case compare(i, 0) of
            GREATER => 1
          | EQUAL => 0
          | LESS => ~1
        infix ==
        infix !=
        fun a == b = a = b
        fun a != b = (a <> b)
        fun *+(b,c,a) = b * c + a
        fun *-(b,c,a) = b * c - b
        fun scan getc = Int.scan StringCvt.DEC getc
    end
structure RNumber : REAL_NUMBER =
    struct
        open Real
        open Real.Math
        type t = Real.real
        val zero = 0.0
        val one = 1.0
        fun signum x = case compare(x,0.0) of
            LESS => ~1.0
          | GREATER => 1.0
          | EQUAL => 0.0
        fun recip x = 1.0 / x
        infix **
        fun i ** n =
            let fun loop 0 = one
                  | loop 1 = i
                  | loop n =
                let val x = loop (Int.div(n, 2))
                    val m = Int.mod(n, 2)
                in
                    if m = 0 then
                        x * x
                    else
                        x * x * i
                end
            in if Int.<(n, 0)
               then raise Domain
               else loop n
            end
        fun max (a, b) = if a < b then b else a
        fun min (a, b) = if a < b then a else b
        fun asinh x = ln (x + sqrt(1.0 + x * x))
        fun acosh x = ln (x + (x + 1.0) * sqrt((x - 1.0)/(x + 1.0)))
        fun atanh x = ln ((1.0 + x) / sqrt(1.0 - x * x))
    end
(*
 Complex(R)     - Functor -
 Provides support for complex numbers based on tuples. Should be
 highly efficient as most operations can be inlined.
 *)
structure CNumber : COMPLEX_NUMBER =
struct
        structure Real = RNumber
        type t = Real.t * Real.t
        type real = Real.t
        val zero = (0.0,0.0)
        val one = (1.0,0.0)
        val pi = (Real.pi, 0.0)
        val e = (Real.e, 0.0)
        fun make (r,i) = (r,i) : t
        fun split z = z
        fun realPart (r,_) = r
        fun imagPart (_,i) = i
        fun abs2 (r,i) = Real.+(Real.*(r,r),Real.*(i,i)) (* FIXME!!! *)
        fun arg (r,i) = Real.atan2(i,r)
        fun modulus z = Real.sqrt(abs2 z)
        fun abs z = (modulus z, 0.0)
        fun signum (z as (r,i)) =
            let val m = modulus z
            in (Real./(r,m), Real./(i,m))
            end
        fun ~ (r1,i1) = (Real.~ r1, Real.~ i1)
        fun (r1,i1) + (r2,i2) = (Real.+(r1,r2), Real.+(i1,i2))
        fun (r1,i1) - (r2,i2) = (Real.-(r1,r2), Real.-(i1,i1))
        fun (r1,i1) * (r2,i2) = (Real.-(Real.*(r1,r2),Real.*(i1,i2)),
                                 Real.+(Real.*(r1,i2),Real.*(r2,i1)))
        fun (r1,i1) / (r2,i2) =
            let val modulus = abs2(r2,i2)
                val (nr,ni) = (r1,i1) * (r2,i2)
            in
                (Real./(nr,modulus), Real./(ni,modulus))
            end
        fun *+((r1,i1),(r2,i2),(r0,i0)) =
            (Real.*+(Real.~ i1, i2, Real.*+(r1,r2,r0)),
             Real.*+(r2, i2, Real.*+(r1,i2,i0)))
        fun *-((r1,i1),(r2,i2),(r0,i0)) =
            (Real.*+(Real.~ i1, i2, Real.*-(r1,r2,r0)),
             Real.*+(r2, i2, Real.*-(r1,i2,i0)))
        infix **
        fun i ** n =
            let fun loop 0 = one
                  | loop 1 = i
                  | loop n =
                let val x = loop (Int.div(n, 2))
                    val m = Int.mod(n, 2)
                in
                    if m = 0 then
                        x * x
                    else
                        x * x * i
                end
            in if Int.<(n, 0)
                   then raise Domain
               else loop n
            end
        fun recip (r1, i1) = 
            let val modulus = abs2(r1, i1)
            in (Real./(r1, modulus), Real./(Real.~ i1, modulus))
            end
        fun ==(z, w) = Real.==(realPart z, realPart w) andalso Real.==(imagPart z, imagPart w)
        fun !=(z, w) = Real.!=(realPart z, realPart w) andalso Real.!=(imagPart z, imagPart w)
        fun fromInt i = (Real.fromInt i, 0.0)
        fun toString (r,i) =
            String.concat ["(",Real.toString r,",",Real.toString i,")"]
        fun exp (x, y) =
            let val expx = Real.exp x
            in (Real.*(x, (Real.cos y)), Real.*(x, (Real.sin y)))
            end
    local
        val half = Real.recip (Real.fromInt 2)
    in
        fun sqrt (z as (x,y)) =
            if Real.==(x, 0.0) andalso Real.==(y, 0.0) then
                zero
            else
                let val m = Real.+(modulus z, Real.abs x)
                    val u' = Real.sqrt (Real.*(m, half))
                    val v' = Real./(Real.abs y , Real.+(u',u'))
                    val (u,v) = if Real.<(x, 0.0) then (v',u') else (u',v')
                in (u, if Real.<(y, 0.0) then Real.~ v else v)
                end
    end
        fun ln z = (Real.ln (modulus z), arg z)
        fun pow (z, n) =
            let val l = ln z
            in exp (l * n)
            end
        fun sin (x, y) = (Real.*(Real.sin x, Real.cosh y),
                          Real.*(Real.cos x, Real.sinh y))
        fun cos (x, y) = (Real.*(Real.cos x, Real.cosh y),
                          Real.~ (Real.*(Real.sin x, Real.sinh y)))
        fun tan (x, y) =
            let val (sx, cx) = (Real.sin x, Real.cos x)
                val (shy, chy) = (Real.sinh y, Real.cosh y)
                val a = (Real.*(sx, chy), Real.*(cx, shy))
                val b = (Real.*(cx, chy), Real.*(Real.~ sx, shy))
            in a / b
            end
        fun sinh (x, y) = (Real.*(Real.cos y, Real.sinh x),
                           Real.*(Real.sin y, Real.cosh x))
        fun cosh (x, y) = (Real.*(Real.cos y, Real.cosh x),
                           Real.*(Real.sin y, Real.sinh x))
        fun tanh (x, y) =
            let val (sy, cy) = (Real.sin y, Real.cos y)
                val (shx, chx) = (Real.sinh x, Real.cosh x)
                val a = (Real.*(cy, shx), Real.*(sy, chx))
                val b = (Real.*(cy, chx), Real.*(sy, shx))
            in a / b
            end
        fun asin (z as (x,y)) =
            let val w = sqrt (one - z * z)
                val (x',y') = ln ((Real.~ y, x) + w)
            in (y', Real.~ x')
            end
        fun acos (z as (x,y)) = 
            let val (x', y') = sqrt (one + z * z)
                val (x'', y'') = ln (z + (Real.~ y', x'))
            in (y'', Real.~ x'')
            end
        fun atan (z as (x,y)) =
            let val w = sqrt (one + z*z)
                val (x',y') = ln ((Real.-(1.0, y), x) / w)
            in (y', Real.~ x')
            end
        fun atan2 (y, x) = atan(y / x)
        fun asinh x = ln (x + sqrt(one + x * x))
        fun acosh x = ln (x + (x + one) * sqrt((x - one)/(x + one)))
        fun atanh x = ln ((one + x) / sqrt(one - x * x))
        fun scan getc =
            let val scanner = Real.scan getc
            in fn stream => 
                  case scanner stream of
                      NONE => NONE
                    | SOME (a, rest) =>
                      case scanner rest of
                          NONE => NONE
                        | SOME (b, rest) => SOME (make(a,b), rest)
            end
end (* ComplexNumber *)
(*
 Copyright (c) Juan Jose Garcia Ripoll.
 All rights reserved.
 Refer to the COPYRIGHT file for license conditions
*)
structure PrettyPrint :>
    sig
        datatype modifier =
                 Int of int |
                 Real of real |
                 Complex of CNumber.t |
                 String of string
        val list : ('a -> string) -> 'a list  -> unit
        val intList : int list -> unit
        val realList : real list -> unit
        val stringList : string list -> unit
        val array : ('a -> string) -> 'a array -> unit
        val intArray : int array -> unit
        val realArray : real array -> unit
        val stringArray : string array -> unit
        val sequence :
            int -> ((int * 'a -> unit) -> 'b -> unit) -> ('a -> string) -> 'b -> unit
        val print : modifier list -> unit
    end =
struct
    datatype modifier =
             Int of int |
             Real of real |
             Complex of CNumber.t |
             String of string     
    fun list _ [] = print "[]"
      | list cvt (a::resta) =
        let fun loop a [] = (print(cvt a); print "]")
              | loop a (b::restb) = (print(cvt a); print ", "; loop b restb)
        in
            print "[";
            loop a resta
        end
    fun boolList a = list Bool.toString a
    fun intList a = list Int.toString a
    fun realList a = list Real.toString a
    fun stringList a = list (fn x => x) a
    fun array cvt a =
        let val length = Array.length a - 1
            fun print_one (i,x) =
                (print(cvt x); if not(i = length) then print ", " else ())
        in
            Array.appi print_one a
        end
    fun boolArray a = array Bool.toString a
    fun intArray a = array Int.toString a
    fun realArray a = array Real.toString a
    fun stringArray a = array (fn x => x) a
    fun sequence length appi cvt seq =
        let val length = length - 1
            fun print_one (i:int,x) =
                (print(cvt x); if not(i = length) then print ", " else ())
        in
            print "[";
            appi print_one seq;
            print "]\n"
        end
    fun print b =
        let fun printer (Int a) = INumber.toString a
              | printer (Real a) = RNumber.toString a
              | printer (Complex a) = CNumber.toString a
              | printer (String a) = a
        in  List.app (fn x => (TextIO.print (printer x))) b
        end
end (* PrettyPrint *)
fun print' x = List.app print x
(*
 Copyright (c) Juan Jose Garcia Ripoll.
 All rights reserved.
 Refer to the COPYRIGHT file for license conditions
*)
structure INumberArray =
    struct
        open Array
        type array = INumber.t array
        type vector = INumber.t vector
        type elem  = INumber.t
        structure Vector =
            struct
                open Vector
                type vector = INumber.t Vector.vector
                type elem = INumber.t
            end
        fun map f a = tabulate(length a, fn x => (f (sub(a,x))))
        fun mapi f a = tabulate(length a, fn x => (f (x,sub(a,x))))
        fun map2 f a b = tabulate(length a, fn x => (f(sub(a,x),sub(b,x))))
    end
structure RNumberArray =
    struct
        open Real64Array
        val sub = Unsafe.Real64Array.sub
        val update = Unsafe.Real64Array.update
        fun map f a = tabulate(length a, fn x => (f (sub(a,x))))
        fun mapi f a = tabulate(length a, fn x => (f (x,sub(a,x))))
        fun map2 f a b = tabulate(length a, fn x => (f(sub(a,x),sub(b,x))))
    end
(*--------------------- COMPLEX ARRAY -------------------------*)
structure BasicCNumberArray =
struct
        structure Complex : COMPLEX_NUMBER = CNumber
        structure Array : MONO_ARRAY = RNumberArray
        type elem = Complex.t
        type array = Array.array * Array.array
        val maxLen = Array.maxLen
        fun length (a,b) = Array.length a
        fun sub ((a,b),index) = Complex.make(Array.sub(a,index),Array.sub(b,index))
        fun update ((a,b),index,z) =
            let val (re,im) = Complex.split z in
                Array.update(a, index, re);
                Array.update(b, index, im)
            end
    local
        fun makeRange (a, start, NONE) = makeRange(a, start, SOME (length a - 1))
          | makeRange (a, start, SOME last) =
            let val len = length a
                val diff = last - start
            in
                if (start >= len) orelse (last >= len) then
                    raise Subscript
                else if diff < 0 then
                    (a, start, 0)
                else
                    (a, start, diff + 1)
            end
    in
        fun array (size,z:elem) =
            let val realsize = size * 2
                val r = Complex.realPart z
                val i = Complex.imagPart z in
                    (Array.array(size,r), Array.array(size,i))
            end
        fun zeroarray size =
            (Array.array(size,Complex.Real.zero),
             Array.array(size,Complex.Real.zero))
        fun tabulate (size,f) =
            let val a = array(size, Complex.zero)
                fun loop i =
                    case i = size of
                        true => a
                      | false => (update(a, i, f i); loop (i+1))
            in
                loop 0
            end
        fun fromList list =
            let val length = List.length list
                val a = zeroarray length
                fun loop (_, []) = a
                  | loop (i, z::rest) = (update(a, i, z);
                                         loop (i+1, rest))
            in
                loop(0,list)
            end
        fun extract range =
            let val (a, start, len) = makeRange range
                fun copy i = sub(a, i + start)
            in tabulate(len, copy)
            end
        fun concat array_list =
            let val total_length = foldl (op +) 0 (map length array_list)
                val a = array(total_length, Complex.zero)
                fun copy (_, []) = a
                  | copy (pos, v::rest) =
                    let fun loop i =
                        case i = 0 of
                            true => ()
                          | false => (update(a, i+pos, sub(v, i)); loop (i-1))
                    in (loop (length v - 1); copy(length v + pos, rest))
                    end
            in
                copy(0, array_list)
            end
        fun copy {src : array, si : int, len : int option, dst : array, di : int } =
            let val (a, ia, la) = makeRange (src, si, len)
                val (b, ib, lb) = makeRange (dst, di, len)
                fun copy i =
                    case i < 0 of
                        true => ()
                      | false => (update(b, i+ib, sub(a, i+ia)); copy (i-1))
            in copy (la - 1)
            end
        val copyVec = copy
        fun modifyi f range =
            let val (a, start, len) = makeRange range
                val last = start + len
                fun loop i =
                    case i >= last of
                        true => ()
                      | false => (update(a, i, f(i, sub(a,i))); loop (i+1))
            in loop start
            end
        fun modify f a =
            let val last = length a
                fun loop i =
                    case i >= last of
                        true => ()
                      | false => (update(a, i, f(sub(a,i))); loop (i+1))
            in loop 0
            end
        fun app f a =
            let val size = length a
                fun loop i =
                    case i = size of
                        true => ()
                      | false => (f(sub(a,i)); loop (i+1))
            in
                loop 0
            end
        fun appi f range =
            let val (a, start, len) = makeRange range
                val last = start + len
                fun loop i =
                    case i >= last of
                        true => ()
                      | false => (f(i, sub(a,i)); loop (i+1))
            in
                loop start
            end
        fun map f a =
            let val len = length a
                val c = zeroarray len
                fun loop ~1 = c
                  | loop i = (update(a, i, f(sub(a,i))); loop (i-1))
            in loop (len-1)
            end
        fun map2 f a b =
            let val len = length a
                val c = zeroarray len
                fun loop ~1 = c
                  | loop i = (update(c, i, f(sub(a,i),sub(b,i)));
                              loop (i-1))
            in loop (len-1)
            end
        fun mapi f range =
            let val (a, start, len) = makeRange range
                fun rule i = f (i+start, sub(a, i+start))
            in tabulate(len, rule)
            end
        fun foldli f init range =
            let val (a, start, len) = makeRange range
                val last = start + len - 1
                fun loop (i, accum) =
                    case i > last of
                        true => accum
                      | false => loop (i+1, f(i, sub(a,i), accum))
            in loop (start, init)
            end
        fun foldri f init range =
            let val (a, start, len) = makeRange range
                val last = start + len - 1
                fun loop (i, accum) =
                    case i < start of
                        true => accum
                      | false => loop (i-1, f(i, sub(a,i), accum))
            in loop (last, init)
            end
        fun foldl f init a = foldli (fn (_, a, x) => f(a,x)) init (a,0,NONE)
        fun foldr f init a = foldri (fn (_, x, a) => f(x,a)) init (a,0,NONE)
    end
end (* BasicCNumberArray *)
structure CNumberArray =
    struct
        structure Vector =
            struct
                open BasicCNumberArray
                type vector = array
            end : MONO_VECTOR
        type vector = Vector.vector
        open BasicCNumberArray
    end (* CNumberArray *)
structure ITensor =
    struct
        structure Number = INumber
        structure Array = INumberArray
(*
 Copyright (c) Juan Jose Garcia Ripoll.
 All rights reserved.
 Refer to the COPYRIGHT file for license conditions
*)
structure MonoTensor  =
    struct
(* PARAMETERS
        structure Array = Array
*)
        structure Index  = Index
        type elem = Array.elem
        type index = Index.t
        type tensor = {shape : index, indexer : Index.indexer, data : Array.array}
        type t = tensor
        exception Shape
        exception Match
        exception Index
    local
    (*----- LOCALS -----*)
        fun make' (shape, data) =
            {shape = shape, indexer = Index.indexer shape, data = data}
        fun toInt {shape, indexer, data} index = indexer index
        fun splitList (l as (a::rest), place) =
            let fun loop (left,here,right) 0 =  (List.rev left,here,right)
                  | loop (_,_,[]) place = raise Index
                  | loop (left,here,a::right) place = 
                loop (here::left,a,right) (place-1)
            in
                if place <= 0 then
                    loop ([],a,rest) (List.length rest - place)
                else
                    loop ([],a,rest) (place - 1)
            end
    in
    (*----- STRUCTURAL OPERATIONS & QUERIES ------*)
        fun new (shape, init) =
            if not (Index.validShape shape) then
                raise Shape
            else
                let val length = Index.length shape in
                    {shape = shape,
                     indexer = Index.indexer shape,
                     data = Array.array(length,init)}
                end
        fun toArray {shape, indexer, data} = data
        fun length {shape, indexer, data} =  Array.length data
        fun shape {shape, indexer, data} = shape
        fun rank t = List.length (shape t)
        fun reshape new_shape tensor =
            if Index.validShape new_shape then
                case (Index.length new_shape) = length tensor of
                    true => make'(new_shape, toArray tensor)
                  | false => raise Match
            else
                raise Shape
        fun fromArray (s, a) =
            case Index.validShape s andalso 
                 ((Index.length s) = (Array.length a)) of
                 true => make'(s, a)
               | false => raise Shape
        fun fromList (s, a) = fromArray (s, Array.fromList a)
        fun tabulate (shape,f) =
            if Index.validShape shape then
                let val last = Index.last shape
                    val length = Index.length shape
                    val c = Array.array(length, f last)
                    fun dotable (c, indices, i) =
                        (Array.update(c, i, f indices);
                         if i <= 1
                         then c
                         else dotable(c, Index.prev' shape indices, i-1))
                in make'(shape,dotable(c, Index.prev' shape last, length-2))
                end
            else
                raise Shape
        (*----- ELEMENTWISE OPERATIONS -----*)
        fun sub (t, index) = Array.sub(#data t, toInt t index)
        fun update (t, index, value) =
            Array.update(toArray t, toInt t index, value)
        fun map f {shape, indexer, data} =
            {shape = shape, indexer = indexer, data = Array.map f data}
        fun map2 f t1 t2=
            let val {shape=shape1, indexer=indexer1, data=data1} = t1
                val {shape=shape2, indexer=indexer2, data=data2} = t2
            in
                if Index.eq(shape1,shape2) then
                    {shape = shape1,
                     indexer = indexer1,
                     data = Array.map2 f data1 data2}
                else
                    raise Match
        end
        fun appi f tensor = Array.appi f (toArray tensor)
        fun app f tensor = Array.app f (toArray tensor)
        fun all f tensor =
            let val a = toArray tensor
            in Loop.all(0, length tensor - 1, fn i =>
                        f (Array.sub(a, i)))
            end
        fun any f tensor =
            let val a = toArray tensor
            in Loop.any(0, length tensor - 1, fn i =>
                        f (Array.sub(a, i)))
            end
        fun foldl f init tensor = Array.foldl f init (toArray tensor)
        fun foldln f init {shape, indexer, data=a} index =
            let val (head,lk,tail) = splitList(shape, index)
                val li = Index.length head
                val lj = Index.length tail
                val c = Array.array(li * lj,init)
                fun loopi (0, _,  _)  = ()
                  | loopi (i, ia, ic) =
                    (Array.update(c, ic, f(Array.sub(c,ic), Array.sub(a,ia)));
                     loopi (i-1, ia+1, ic+1))
                fun loopk (0, ia, _)  = ia
                  | loopk (k, ia, ic) = (loopi (li, ia, ic);
                                         loopk (k-1, ia+li, ic))
                fun loopj (0, _,  _)  = ()
                  | loopj (j, ia, ic) = loopj (j-1, loopk(lk,ia,ic), ic+li)
            in
                loopj (lj, 0, 0);
                make'(head @ tail, c)
            end
        (* --- POLYMORPHIC ELEMENTWISE OPERATIONS --- *)
        fun array_map' f a =
            let fun apply index = f(Array.sub(a,index)) in
                Tensor.Array.tabulate(Array.length a, apply)
            end
        fun map' f t = Tensor.fromArray(shape t, array_map' f (toArray t))
        fun map2' f t1 t2 =
            let val d1 = toArray t1
                val d2 = toArray t2
                fun apply i = f (Array.sub(d1,i), Array.sub(d2,i))
                val len = Array.length d1
            in
                if Index.eq(shape t1, shape t2) then
                    Tensor.fromArray(shape t1, Tensor.Array.tabulate(len,apply))
                else
                    raise Match
            end
        fun foldl' f init {shape, indexer, data=a} index =
            let val (head,lk,tail) = splitList(shape, index)
                val li = Index.length head
                val lj = Index.length tail
                val c = Tensor.Array.array(li * lj,init)
                fun loopi (0, _,  _)  = ()
                  | loopi (i, ia, ic) =
                    (Tensor.Array.update(c,ic,f(Tensor.Array.sub(c,ic),Array.sub(a,ia)));
                     loopi (i-1, ia+1, ic+1))
                fun loopk (0, ia, _)  = ia
                  | loopk (k, ia, ic) = (loopi (li, ia, ic);
                                         loopk (k-1, ia+li, ic))
                fun loopj (0, _,  _)  = ()
                  | loopj (j, ia, ic) = loopj (j-1, loopk(lk,ia,ic), ic+li)
            in
                loopj (lj, 0, 0);
                make'(head @ tail, c)
            end
    end
    end (* MonoTensor *)
        open MonoTensor
    local
        (*
         LEFT INDEX CONTRACTION:
         a = a(i1,i2,...,in)
         b = b(j1,j2,...,jn)
         c = c(i2,...,in,j2,...,jn)
         = sum(a(k,i2,...,jn)*b(k,j2,...jn)) forall k
         MEANINGFUL VARIABLES:
         lk = i1 = j1
         li = i2*...*in
         lj = j2*...*jn
         *)
        fun do_fold_first a b c lk lj li =
            let fun loopk (0, _,  _,  accum) = accum
                  | loopk (k, ia, ib, accum) =
                    let val delta = Number.*(Array.sub(a,ia),Array.sub(b,ib))
                    in loopk (k-1, ia+1, ib+1, Number.+(delta,accum))
                    end
                fun loopj (0, ib, ic) = c
                  | loopj (j, ib, ic) =
                    let fun loopi (0, ia, ic) = ic
                          | loopi (i, ia, ic) =
                        (Array.update(c, ic, loopk(lk, ia, ib, Number.zero));
                         loopi(i-1, ia+lk, ic+1))
                    in
                        loopj(j-1, ib+lk, loopi(li, 0, ic))
                    end
            in loopj(lj, 0, 0)
            end
    in
        fun +* ta tb =
            let val (rank_a,lk::rest_a,a) = (rank ta, shape ta, toArray ta)
                val (rank_b,lk2::rest_b,b) = (rank tb, shape tb, toArray tb)
            in if not(lk = lk2)
               then raise Match
               else let val li = Index.length rest_a
                        val lj = Index.length rest_b
                        val c = Array.array(li*lj,Number.zero)
                    in fromArray(rest_a @ rest_b,
                                 do_fold_first a b c lk li lj)
                    end
            end
    end
    local
        (*
         LAST INDEX CONTRACTION:
         a = a(i1,i2,...,in)
         b = b(j1,j2,...,jn)
         c = c(i2,...,in,j2,...,jn)
         = sum(mult(a(i1,i2,...,k),b(j1,j2,...,k))) forall k
         MEANINGFUL VARIABLES:
         lk = in = jn
         li = i1*...*i(n-1)
         lj = j1*...*j(n-1)
         *)
        fun do_fold_last a b c lk lj li =
            let fun loopi (0, ia, ic, fac) = ()
                  | loopi (i, ia, ic, fac) =
                    let val old = Array.sub(c,ic)
                        val inc = Number.*(Array.sub(a,ia),fac)
                    in
                        Array.update(c,ic,Number.+(old,inc));
                        loopi(i-1, ia+1, ic+1, fac)
                    end
                fun loopj (j, ib, ic) =
                    let fun loopk (0, ia, ib) = ()
                          | loopk (k, ia, ib) =
                            (loopi(li, ia, ic, Array.sub(b,ib));
                             loopk(k-1, ia+li, ib+lj))
                    in case j of
                           0 => c
                         | _ => (loopk(lk, 0, ib);
                                 loopj(j-1, ib+1, ic+li))
                    end (* loopj *)
            in
                loopj(lj, 0, 0)
            end
    in
        fun *+ ta tb  =
            let val (rank_a,shape_a,a) = (rank ta, shape ta, toArray ta)
                val (rank_b,shape_b,b) = (rank tb, shape tb, toArray tb)
                val (lk::rest_a) = List.rev shape_a
                val (lk2::rest_b) = List.rev shape_b
            in if not(lk = lk2)
               then raise Match
               else let val li = Index.length rest_a
                        val lj = Index.length rest_b
                        val c = Array.array(li*lj,Number.zero)
                    in fromArray(List.rev rest_a @ List.rev rest_b,
                                 do_fold_last a b c lk li lj)
                    end
            end
    end
        (* ALGEBRAIC OPERATIONS *)
        infix **
        infix ==
        infix !=
        fun a + b = map2 Number.+ a b
        fun a - b = map2 Number.- a b
        fun a * b = map2 Number.* a b
        fun a ** i = map (fn x => (Number.**(x,i))) a
        fun ~ a = map Number.~ a
        fun abs a = map Number.abs a
        fun signum a = map Number.signum a
        fun a == b = map2' Number.== a b
        fun a != b = map2' Number.!= a b
        fun toString a = raise Domain
        fun fromInt a = new([1], Number.fromInt a)
        (* TENSOR SPECIFIC OPERATIONS *)
        fun *> n = map (fn x => Number.*(n,x))
        fun print t =
            (PrettyPrint.intList (shape t);
             TextIO.print "\n";
             PrettyPrint.sequence (length t) appi Number.toString t)
        fun normInf a =
            let fun accum (y,x) = Number.max(x,Number.abs y)
            in  foldl accum Number.zero a
            end
    end (* NumberTensor *)
structure RTensor =
    struct
        structure Number = RNumber
        structure Array = RNumberArray
(*
 Copyright (c) Juan Jose Garcia Ripoll.
 All rights reserved.
 Refer to the COPYRIGHT file for license conditions
*)
structure MonoTensor  =
    struct
(* PARAMETERS
        structure Array = Array
*)
        structure Index  = Index
        type elem = Array.elem
        type index = Index.t
        type tensor = {shape : index, indexer : Index.indexer, data : Array.array}
        type t = tensor
        exception Shape
        exception Match
        exception Index
    local
    (*----- LOCALS -----*)
        fun make' (shape, data) =
            {shape = shape, indexer = Index.indexer shape, data = data}
        fun toInt {shape, indexer, data} index = indexer index
        fun splitList (l as (a::rest), place) =
            let fun loop (left,here,right) 0 =  (List.rev left,here,right)
                  | loop (_,_,[]) place = raise Index
                  | loop (left,here,a::right) place = 
                loop (here::left,a,right) (place-1)
            in
                if place <= 0 then
                    loop ([],a,rest) (List.length rest - place)
                else
                    loop ([],a,rest) (place - 1)
            end
    in
    (*----- STRUCTURAL OPERATIONS & QUERIES ------*)
        fun new (shape, init) =
            if not (Index.validShape shape) then
                raise Shape
            else
                let val length = Index.length shape in
                    {shape = shape,
                     indexer = Index.indexer shape,
                     data = Array.array(length,init)}
                end
        fun toArray {shape, indexer, data} = data
        fun length {shape, indexer, data} =  Array.length data
        fun shape {shape, indexer, data} = shape
        fun rank t = List.length (shape t)
        fun reshape new_shape tensor =
            if Index.validShape new_shape then
                case (Index.length new_shape) = length tensor of
                    true => make'(new_shape, toArray tensor)
                  | false => raise Match
            else
                raise Shape
        fun fromArray (s, a) =
            case Index.validShape s andalso 
                 ((Index.length s) = (Array.length a)) of
                 true => make'(s, a)
               | false => raise Shape
        fun fromList (s, a) = fromArray (s, Array.fromList a)
        fun tabulate (shape,f) =
            if Index.validShape shape then
                let val last = Index.last shape
                    val length = Index.length shape
                    val c = Array.array(length, f last)
                    fun dotable (c, indices, i) =
                        (Array.update(c, i, f indices);
                         if i <= 1
                         then c
                         else dotable(c, Index.prev' shape indices, i-1))
                in make'(shape,dotable(c, Index.prev' shape last, length-2))
                end
            else
                raise Shape
        (*----- ELEMENTWISE OPERATIONS -----*)
        fun sub (t, index) = Array.sub(#data t, toInt t index)
        fun update (t, index, value) =
            Array.update(toArray t, toInt t index, value)
        fun map f {shape, indexer, data} =
            {shape = shape, indexer = indexer, data = Array.map f data}
        fun map2 f t1 t2=
            let val {shape=shape1, indexer=indexer1, data=data1} = t1
                val {shape=shape2, indexer=indexer2, data=data2} = t2
            in
                if Index.eq(shape1,shape2) then
                    {shape = shape1,
                     indexer = indexer1,
                     data = Array.map2 f data1 data2}
                else
                    raise Match
        end
        fun appi f tensor = Array.appi f (toArray tensor)
        fun app f tensor = Array.app f (toArray tensor)
        fun all f tensor =
            let val a = toArray tensor
            in Loop.all(0, length tensor - 1, fn i =>
                        f (Array.sub(a, i)))
            end
        fun any f tensor =
            let val a = toArray tensor
            in Loop.any(0, length tensor - 1, fn i =>
                        f (Array.sub(a, i)))
            end
        fun foldl f init tensor = Array.foldl f init (toArray tensor)
        fun foldln f init {shape, indexer, data=a} index =
            let val (head,lk,tail) = splitList(shape, index)
                val li = Index.length head
                val lj = Index.length tail
                val c = Array.array(li * lj,init)
                fun loopi (0, _,  _)  = ()
                  | loopi (i, ia, ic) =
                    (Array.update(c, ic, f(Array.sub(c,ic), Array.sub(a,ia)));
                     loopi (i-1, ia+1, ic+1))
                fun loopk (0, ia, _)  = ia
                  | loopk (k, ia, ic) = (loopi (li, ia, ic);
                                         loopk (k-1, ia+li, ic))
                fun loopj (0, _,  _)  = ()
                  | loopj (j, ia, ic) = loopj (j-1, loopk(lk,ia,ic), ic+li)
            in
                loopj (lj, 0, 0);
                make'(head @ tail, c)
            end
        (* --- POLYMORPHIC ELEMENTWISE OPERATIONS --- *)
        fun array_map' f a =
            let fun apply index = f(Array.sub(a,index)) in
                Tensor.Array.tabulate(Array.length a, apply)
            end
        fun map' f t = Tensor.fromArray(shape t, array_map' f (toArray t))
        fun map2' f t1 t2 =
            let val d1 = toArray t1
                val d2 = toArray t2
                fun apply i = f (Array.sub(d1,i), Array.sub(d2,i))
                val len = Array.length d1
            in
                if Index.eq(shape t1, shape t2) then
                    Tensor.fromArray(shape t1, Tensor.Array.tabulate(len,apply))
                else
                    raise Match
            end
        fun foldl' f init {shape, indexer, data=a} index =
            let val (head,lk,tail) = splitList(shape, index)
                val li = Index.length head
                val lj = Index.length tail
                val c = Tensor.Array.array(li * lj,init)
                fun loopi (0, _,  _)  = ()
                  | loopi (i, ia, ic) =
                    (Tensor.Array.update(c,ic,f(Tensor.Array.sub(c,ic),Array.sub(a,ia)));
                     loopi (i-1, ia+1, ic+1))
                fun loopk (0, ia, _)  = ia
                  | loopk (k, ia, ic) = (loopi (li, ia, ic);
                                         loopk (k-1, ia+li, ic))
                fun loopj (0, _,  _)  = ()
                  | loopj (j, ia, ic) = loopj (j-1, loopk(lk,ia,ic), ic+li)
            in
                loopj (lj, 0, 0);
                make'(head @ tail, c)
            end
    end
    end (* MonoTensor *)
        open MonoTensor
    local
        (*
         LEFT INDEX CONTRACTION:
         a = a(i1,i2,...,in)
         b = b(j1,j2,...,jn)
         c = c(i2,...,in,j2,...,jn)
         = sum(a(k,i2,...,jn)*b(k,j2,...jn)) forall k
         MEANINGFUL VARIABLES:
         lk = i1 = j1
         li = i2*...*in
         lj = j2*...*jn
         *)
        fun do_fold_first a b c lk lj li =
            let fun loopk (0, _,  _,  accum) = accum
                  | loopk (k, ia, ib, accum) =
                    let val delta = Number.*(Array.sub(a,ia),Array.sub(b,ib))
                    in loopk (k-1, ia+1, ib+1, Number.+(delta,accum))
                    end
                fun loopj (0, ib, ic) = c
                  | loopj (j, ib, ic) =
                    let fun loopi (0, ia, ic) = ic
                          | loopi (i, ia, ic) =
                        (Array.update(c, ic, loopk(lk, ia, ib, Number.zero));
                         loopi(i-1, ia+lk, ic+1))
                    in
                        loopj(j-1, ib+lk, loopi(li, 0, ic))
                    end
            in loopj(lj, 0, 0)
            end
    in
        fun +* ta tb =
            let val (rank_a,lk::rest_a,a) = (rank ta, shape ta, toArray ta)
                val (rank_b,lk2::rest_b,b) = (rank tb, shape tb, toArray tb)
            in if not(lk = lk2)
               then raise Match
               else let val li = Index.length rest_a
                        val lj = Index.length rest_b
                        val c = Array.array(li*lj,Number.zero)
                    in fromArray(rest_a @ rest_b,
                                 do_fold_first a b c lk li lj)
                    end
            end
    end
    local
        (*
         LAST INDEX CONTRACTION:
         a = a(i1,i2,...,in)
         b = b(j1,j2,...,jn)
         c = c(i2,...,in,j2,...,jn)
         = sum(mult(a(i1,i2,...,k),b(j1,j2,...,k))) forall k
         MEANINGFUL VARIABLES:
         lk = in = jn
         li = i1*...*i(n-1)
         lj = j1*...*j(n-1)
         *)
        fun do_fold_last a b c lk lj li =
            let fun loopi (0, ia, ic, fac) = ()
                  | loopi (i, ia, ic, fac) =
                    let val old = Array.sub(c,ic)
                        val inc = Number.*(Array.sub(a,ia),fac)
                    in
                        Array.update(c,ic,Number.+(old,inc));
                        loopi(i-1, ia+1, ic+1, fac)
                    end
                fun loopj (j, ib, ic) =
                    let fun loopk (0, ia, ib) = ()
                          | loopk (k, ia, ib) =
                            (loopi(li, ia, ic, Array.sub(b,ib));
                             loopk(k-1, ia+li, ib+lj))
                    in case j of
                           0 => c
                         | _ => (loopk(lk, 0, ib);
                                 loopj(j-1, ib+1, ic+li))
                    end (* loopj *)
            in
                loopj(lj, 0, 0)
            end
    in
        fun *+ ta tb  =
            let val (rank_a,shape_a,a) = (rank ta, shape ta, toArray ta)
                val (rank_b,shape_b,b) = (rank tb, shape tb, toArray tb)
                val (lk::rest_a) = List.rev shape_a
                val (lk2::rest_b) = List.rev shape_b
            in if not(lk = lk2)
               then raise Match
               else let val li = Index.length rest_a
                        val lj = Index.length rest_b
                        val c = Array.array(li*lj,Number.zero)
                    in fromArray(List.rev rest_a @ List.rev rest_b,
                                 do_fold_last a b c lk li lj)
                    end
            end
    end
        (* ALGEBRAIC OPERATIONS *)
        infix **
        infix ==
        infix !=
        fun a + b = map2 Number.+ a b
        fun a - b = map2 Number.- a b
        fun a * b = map2 Number.* a b
        fun a ** i = map (fn x => (Number.**(x,i))) a
        fun ~ a = map Number.~ a
        fun abs a = map Number.abs a
        fun signum a = map Number.signum a
        fun a == b = map2' Number.== a b
        fun a != b = map2' Number.!= a b
        fun toString a = raise Domain
        fun fromInt a = new([1], Number.fromInt a)
        (* TENSOR SPECIFIC OPERATIONS *)
        fun *> n = map (fn x => Number.*(n,x))
        fun print t =
            (PrettyPrint.intList (shape t);
             TextIO.print "\n";
             PrettyPrint.sequence (length t) appi Number.toString t)
        fun a / b = map2 Number./ a b
        fun recip a = map Number.recip a
        fun ln a = map Number.ln a
        fun pow (a, b) = map (fn x => (Number.pow(x,b))) a
        fun exp a = map Number.exp a
        fun sqrt a = map Number.sqrt a
        fun cos a = map Number.cos a
        fun sin a = map Number.sin a
        fun tan a = map Number.tan a
        fun sinh a = map Number.sinh a
        fun cosh a = map Number.cosh a
        fun tanh a = map Number.tanh a
        fun asin a = map Number.asin a
        fun acos a = map Number.acos a
        fun atan a = map Number.atan a
        fun asinh a = map Number.asinh a
        fun acosh a = map Number.acosh a
        fun atanh a = map Number.atanh a
        fun atan2 (a,b) = map2 Number.atan2 a b
        fun normInf a =
            let fun accum (y,x) = Number.max(x,Number.abs y)
            in  foldl accum Number.zero a
            end
        fun norm1 a =
            let fun accum (y,x) = Number.+(x,Number.abs y)
            in  foldl accum Number.zero a
            end
        fun norm2 a =
            let fun accum (y,x) = Number.+(x, Number.*(y,y))
            in Number.sqrt(foldl accum Number.zero a)
            end
    end (* RTensor *)
structure CTensor =
struct
    structure Number = CNumber
    structure Array = CNumberArray
(*
 Copyright (c) Juan Jose Garcia Ripoll.
 All rights reserved.
 Refer to the COPYRIGHT file for license conditions
*)
structure MonoTensor  =
    struct
(* PARAMETERS
        structure Array = Array
*)
        structure Index  = Index
        type elem = Array.elem
        type index = Index.t
        type tensor = {shape : index, indexer : Index.indexer, data : Array.array}
        type t = tensor
        exception Shape
        exception Match
        exception Index
    local
    (*----- LOCALS -----*)
        fun make' (shape, data) =
            {shape = shape, indexer = Index.indexer shape, data = data}
        fun toInt {shape, indexer, data} index = indexer index
        fun splitList (l as (a::rest), place) =
            let fun loop (left,here,right) 0 =  (List.rev left,here,right)
                  | loop (_,_,[]) place = raise Index
                  | loop (left,here,a::right) place = 
                loop (here::left,a,right) (place-1)
            in
                if place <= 0 then
                    loop ([],a,rest) (List.length rest - place)
                else
                    loop ([],a,rest) (place - 1)
            end
    in
    (*----- STRUCTURAL OPERATIONS & QUERIES ------*)
        fun new (shape, init) =
            if not (Index.validShape shape) then
                raise Shape
            else
                let val length = Index.length shape in
                    {shape = shape,
                     indexer = Index.indexer shape,
                     data = Array.array(length,init)}
                end
        fun toArray {shape, indexer, data} = data
        fun length {shape, indexer, data} =  Array.length data
        fun shape {shape, indexer, data} = shape
        fun rank t = List.length (shape t)
        fun reshape new_shape tensor =
            if Index.validShape new_shape then
                case (Index.length new_shape) = length tensor of
                    true => make'(new_shape, toArray tensor)
                  | false => raise Match
            else
                raise Shape
        fun fromArray (s, a) =
            case Index.validShape s andalso 
                 ((Index.length s) = (Array.length a)) of
                 true => make'(s, a)
               | false => raise Shape
        fun fromList (s, a) = fromArray (s, Array.fromList a)
        fun tabulate (shape,f) =
            if Index.validShape shape then
                let val last = Index.last shape
                    val length = Index.length shape
                    val c = Array.array(length, f last)
                    fun dotable (c, indices, i) =
                        (Array.update(c, i, f indices);
                         if i <= 1
                         then c
                         else dotable(c, Index.prev' shape indices, i-1))
                in make'(shape,dotable(c, Index.prev' shape last, length-2))
                end
            else
                raise Shape
        (*----- ELEMENTWISE OPERATIONS -----*)
        fun sub (t, index) = Array.sub(#data t, toInt t index)
        fun update (t, index, value) =
            Array.update(toArray t, toInt t index, value)
        fun map f {shape, indexer, data} =
            {shape = shape, indexer = indexer, data = Array.map f data}
        fun map2 f t1 t2=
            let val {shape=shape1, indexer=indexer1, data=data1} = t1
                val {shape=shape2, indexer=indexer2, data=data2} = t2
            in
                if Index.eq(shape1,shape2) then
                    {shape = shape1,
                     indexer = indexer1,
                     data = Array.map2 f data1 data2}
                else
                    raise Match
        end
        fun appi f tensor = Array.appi f (toArray tensor, 0, NONE)
        fun app f tensor = Array.app f (toArray tensor)
        fun all f tensor =
            let val a = toArray tensor
            in Loop.all(0, length tensor - 1, fn i =>
                        f (Array.sub(a, i)))
            end
        fun any f tensor =
            let val a = toArray tensor
            in Loop.any(0, length tensor - 1, fn i =>
                        f (Array.sub(a, i)))
            end
        fun foldl f init tensor = Array.foldl f init (toArray tensor)
        fun foldln f init {shape, indexer, data=a} index =
            let val (head,lk,tail) = splitList(shape, index)
                val li = Index.length head
                val lj = Index.length tail
                val c = Array.array(li * lj,init)
                fun loopi (0, _,  _)  = ()
                  | loopi (i, ia, ic) =
                    (Array.update(c, ic, f(Array.sub(c,ic), Array.sub(a,ia)));
                     loopi (i-1, ia+1, ic+1))
                fun loopk (0, ia, _)  = ia
                  | loopk (k, ia, ic) = (loopi (li, ia, ic);
                                         loopk (k-1, ia+li, ic))
                fun loopj (0, _,  _)  = ()
                  | loopj (j, ia, ic) = loopj (j-1, loopk(lk,ia,ic), ic+li)
            in
                loopj (lj, 0, 0);
                make'(head @ tail, c)
            end
        (* --- POLYMORPHIC ELEMENTWISE OPERATIONS --- *)
        fun array_map' f a =
            let fun apply index = f(Array.sub(a,index)) in
                Tensor.Array.tabulate(Array.length a, apply)
            end
        fun map' f t = Tensor.fromArray(shape t, array_map' f (toArray t))
        fun map2' f t1 t2 =
            let val d1 = toArray t1
                val d2 = toArray t2
                fun apply i = f (Array.sub(d1,i), Array.sub(d2,i))
                val len = Array.length d1
            in
                if Index.eq(shape t1, shape t2) then
                    Tensor.fromArray(shape t1, Tensor.Array.tabulate(len,apply))
                else
                    raise Match
            end
        fun foldl' f init {shape, indexer, data=a} index =
            let val (head,lk,tail) = splitList(shape, index)
                val li = Index.length head
                val lj = Index.length tail
                val c = Tensor.Array.array(li * lj,init)
                fun loopi (0, _,  _)  = ()
                  | loopi (i, ia, ic) =
                    (Tensor.Array.update(c,ic,f(Tensor.Array.sub(c,ic),Array.sub(a,ia)));
                     loopi (i-1, ia+1, ic+1))
                fun loopk (0, ia, _)  = ia
                  | loopk (k, ia, ic) = (loopi (li, ia, ic);
                                         loopk (k-1, ia+li, ic))
                fun loopj (0, _,  _)  = ()
                  | loopj (j, ia, ic) = loopj (j-1, loopk(lk,ia,ic), ic+li)
            in
                loopj (lj, 0, 0);
                make'(head @ tail, c)
            end
    end
    end (* MonoTensor *)
    open MonoTensor
    local
        (*
         LEFT INDEX CONTRACTION:
         a = a(i1,i2,...,in)
         b = b(j1,j2,...,jn)
         c = c(i2,...,in,j2,...,jn)
         = sum(a(k,i2,...,jn)*b(k,j2,...jn)) forall k
         MEANINGFUL VARIABLES:
         lk = i1 = j1
         li = i2*...*in
         lj = j2*...*jn
         *)
        fun do_fold_first a b c lk lj li =
            let fun loopk (0, _, _, r, i) = Number.make(r,i)
                  | loopk (k, ia, ib, r, i) =
                    let val (ar, ai) = Array.sub(a,ia)
                        val (br, bi) = Array.sub(b,ib)
                        val dr = ar * br - ai * bi
                        val di = ar * bi + ai * br
                    in loopk (k-1, ia+1, ib+1, r+dr, i+di)
                    end
                fun loopj (0, ib, ic) = c
                  | loopj (j, ib, ic) =
                    let fun loopi (0, ia, ic) = ic
                          | loopi (i, ia, ic) =
                            (Array.update(c, ic, loopk(lk, ia, ib, RNumber.zero, RNumber.zero));
                             loopi(i-1, ia+lk, ic+1))
                    in loopj(j-1, ib+lk, loopi(li, 0, ic))
                    end
            in loopj(lj, 0, 0)
            end
    in
        fun +* ta tb =
            let val (rank_a,lk::rest_a,a) = (rank ta, shape ta, toArray ta)
                val (rank_b,lk2::rest_b,b) = (rank tb, shape tb, toArray tb)
            in if not(lk = lk2)
               then raise Match
               else let val li = Index.length rest_a
                        val lj = Index.length rest_b
                        val c = Array.array(li*lj,Number.zero)
                    in fromArray(rest_a @ rest_b, do_fold_first a b c lk li lj)
                    end
            end
    end
    local
        (*
         LAST INDEX CONTRACTION:
         a = a(i1,i2,...,in)
         b = b(j1,j2,...,jn)
         c = c(i2,...,in,j2,...,jn)
         = sum(mult(a(i1,i2,...,k),b(j1,j2,...,k))) forall k
         MEANINGFUL VARIABLES:
         lk = in = jn
         li = i1*...*i(n-1)
         lj = j1*...*j(n-1)
         *)
        fun do_fold_last a b c lk lj li =
            let fun loopi(0, _, _, _, _) = ()
                  | loopi(i, ia, ic, br, bi) =
                    let val (cr,ci) = Array.sub(c,ic)
                        val (ar,ai) = Array.sub(a,ia)
                        val dr = (ar * br - ai * bi)
                        val di = (ar * bi + ai * br)
                    in
                        Array.update(c,ic,Number.make(cr+dr,ci+di));
                        loopi(i-1, ia+1, ic+1, br, bi)
                    end
                fun loopj(j, ib, ic) =
                    let fun loopk(0, _, _) = ()
                          | loopk(k, ia, ib) =
                            let val (br, bi) = Array.sub(b,ib)
                            in
                                loopi(li, ia, ic, br, bi);
                                loopk(k-1, ia+li, ib+lj)
                            end
                in case j of
                    0 => c
                  | _ => (loopk(lk, 0, ib);
                          loopj(j-1, ib+1, ic+li))
                end (* loopj *)
            in
                loopj(lj, 0, 0)
            end
    in
        fun *+ ta tb  =
            let val (rank_a,shape_a,a) = (rank ta, shape ta, toArray ta)
                val (rank_b,shape_b,b) = (rank tb, shape tb, toArray tb)
                val (lk::rest_a) = List.rev shape_a
                val (lk2::rest_b) = List.rev shape_b
            in
                if not(lk = lk2) then
                    raise Match
                else
                    let val li = Index.length rest_a
                        val lj = Index.length rest_b
                        val c = Array.array(li*lj,Number.zero)
                    in
                        fromArray(List.rev rest_a @ List.rev rest_b,
                                  do_fold_last a b c lk li lj)
                    end
            end
    end
    (* ALGEBRAIC OPERATIONS *)
    infix **
    infix ==
    infix !=
    fun a + b = map2 Number.+ a b
    fun a - b = map2 Number.- a b
    fun a * b = map2 Number.* a b
    fun a ** i = map (fn x => (Number.**(x,i))) a
    fun ~ a = map Number.~ a
    fun abs a = map Number.abs a
    fun signum a = map Number.signum a
    fun a == b = map2' Number.== a b
    fun a != b = map2' Number.!= a b
    fun toString a = raise Domain
    fun fromInt a = new([1], Number.fromInt a)
    (* TENSOR SPECIFIC OPERATIONS *)
    fun *> n = map (fn x => Number.*(n,x))
    fun print t =
        (PrettyPrint.intList (shape t);
         TextIO.print "\n";
         PrettyPrint.sequence (length t) appi Number.toString t)
    fun a / b = map2 Number./ a b
    fun recip a = map Number.recip a
    fun ln a = map Number.ln a
    fun pow (a, b) = map (fn x => (Number.pow(x,b))) a
    fun exp a = map Number.exp a
    fun sqrt a = map Number.sqrt a
    fun cos a = map Number.cos a
    fun sin a = map Number.sin a
    fun tan a = map Number.tan a
    fun sinh a = map Number.sinh a
    fun cosh a = map Number.cosh a
    fun tanh a = map Number.tanh a
    fun asin a = map Number.asin a
    fun acos a = map Number.acos a
    fun atan a = map Number.atan a
    fun asinh a = map Number.asinh a
    fun acosh a = map Number.acosh a
    fun atanh a = map Number.atanh a
    fun atan2 (a,b) = map2 Number.atan2 a b
    fun normInf a =
        let fun accum (y,x) = RNumber.max(x, Number.realPart(Number.abs y))
        in  foldl accum RNumber.zero a
        end
    fun norm1 a =
        let fun accum (y,x) = RNumber.+(x, Number.realPart(Number.abs y))
        in  foldl accum RNumber.zero a
        end
    fun norm2 a =
        let fun accum (y,x) = RNumber.+(x, Number.abs2 y)
        in RNumber.sqrt(foldl accum RNumber.zero a)
        end
end (* CTensor *)
structure MathFile =
struct

type file = TextIO.instream

exception Data

fun assert NONE = raise Data
  | assert (SOME a) = a

(* ------------------ INPUT --------------------- *)

fun intRead file = assert(TextIO.scanStream INumber.scan file)
fun realRead file = assert(TextIO.scanStream RNumber.scan file)
fun complexRead file = assert(TextIO.scanStream CNumber.scan file)

fun listRead eltScan file =
    let val length = intRead file
        fun eltRead file = assert(TextIO.scanStream eltScan file)
        fun loop (0,accum) = accum
          | loop (i,accum) = loop(i-1, eltRead file :: accum)
    in
        if length < 0
        then raise Data
        else List.rev(loop(length,[]))
    end

fun intListRead file = listRead INumber.scan file
fun realListRead file = listRead RNumber.scan file
fun complexListRead file = listRead CNumber.scan file

fun intTensorRead file =
    let val shape = intListRead file
        val length = Index.length shape
        val first = intRead file
        val a = ITensor.Array.array(length, first)
        fun loop 0 = ITensor.fromArray(shape, a)
          | loop j = (ITensor.Array.update(a, length-j, intRead file);
                      loop (j-1))
    in loop (length - 1)
    end

fun realTensorRead file =
    let val shape = intListRead file
        val length = Index.length shape
        val first = realRead file
        val a = RTensor.Array.array(length, first)
        fun loop 0 = RTensor.fromArray(shape, a)
          | loop j = (RTensor.Array.update(a, length-j, realRead file);
                      loop (j-1))
    in loop (length - 1)
    end

fun complexTensorRead file =
    let val shape = intListRead file
        val length = Index.length shape
        val first = complexRead file
        val a = CTensor.Array.array(length, first)
        fun loop j = if j = length
                     then CTensor.fromArray(shape, a)
                     else (CTensor.Array.update(a, j, complexRead file);
                           loop (j+1))
    in loop 1
    end

(* ------------------ OUTPUT -------------------- *)
fun linedOutput(file, x) = (TextIO.output(file, x); TextIO.output(file, "\n"))

fun intWrite file x = linedOutput(file, INumber.toString x)
fun realWrite file x = linedOutput(file, RNumber.toString x)
fun complexWrite file x =
    let val (r,i) = CNumber.split x
    in linedOutput(file, concat [RNumber.toString r, " ", RNumber.toString i])
    end

fun listWrite converter file x =
    (intWrite file (length x);
     List.app (fn x => (linedOutput(file, converter x))) x)

fun intListWrite file x = listWrite INumber.toString file x
fun realListWrite file x = listWrite RNumber.toString file x
fun complexListWrite file x = listWrite CNumber.toString file x

fun intTensorWrite file x = (intListWrite file (ITensor.shape x); ITensor.app (fn x => (intWrite file x)) x)
fun realTensorWrite file x = (intListWrite file (RTensor.shape x); RTensor.app (fn x => (realWrite file x)) x)
fun complexTensorWrite file x = (intListWrite file (CTensor.shape x); CTensor.app (fn x => (complexWrite file x)) x)
end

fun loop 0 _ = ()
  | loop n f = (f(); loop (n-1) f)

fun test_operator new list_op list_sizes =
    let fun test_many list_op size =
        let fun test_op (times,f) =
            let val a = new size
            in (EvalTimer.timerOn();
                loop times (fn _ => f(a,a));
                let val t = LargeInt.toInt(EvalTimer.timerRead()) div times
                    val i = StringCvt.padLeft #" " 6 (Int.toString t)
                in print i
                end)
            end
        in
            print (Int.toString size);
            print " ";
            List.app test_op list_op;
            print "\n"
        end
    in List.app (test_many list_op) list_sizes
    end

structure Main =
   struct
      fun one() =
         let
            val _ =
               let val operators = [(20, RTensor.+), (20, RTensor.* ), (20, RTensor./),
                                    (4, fn (a,b) => RTensor.+* a b),
                                    (4, fn (a,b) => RTensor.*+ a b)]
                  fun constructor size = RTensor.new([size,size],1.0)
               in
                  print "Real tensors: (+, *, /, +*, *+)\n";
                  test_operator constructor operators [100,200,300,400,500];
                  print "\n\n"
               end
            
      val _ =
         let val operators = [(20, CTensor.+), (20, CTensor.* ), (20, CTensor./),
                              (4, fn (a,b) => CTensor.+* a b),
                              (4, fn (a,b) => CTensor.*+ a b)]
            fun constructor size = CTensor.new([size,size],CNumber.one)
         in
            print "Real tensors: (+, *, /, +*, *+)\n";
            test_operator constructor operators [100,200,300,400,500];
            print "\n\n"
         end
         in ()
         end

      fun doit n =
         if n = 0
            then ()
         else (one ()
               ; doit (n - 1))
   end