1 (* Obtained at http
://www
.arrakis
.es
/~worm
/ *)
3 signature MONO_VECTOR
=
8 val fromList
: elem list
-> vector
9 val tabulate
: (int * (int -> elem
)) -> vector
10 val length
: vector
-> int
11 val sub
: (vector
* int) -> elem
12 val extract
: (vector
* int * int option
) -> vector
13 val concat
: vector list
-> vector
14 val mapi
: ((int * elem
) -> elem
) -> (vector
* int * int option
) -> vector
15 val map
: (elem
-> elem
) -> vector
-> vector
16 val appi
: ((int * elem
) -> unit
) -> (vector
* int * int option
) -> unit
17 val app
: (elem
-> unit
) -> vector
-> unit
18 val foldli
: ((int * elem
* 'a
) -> 'a
) -> 'a
-> (vector
* int * int option
) -> 'a
19 val foldri
: ((int * elem
* 'a
) -> 'a
) -> 'a
-> (vector
* int * int option
) -> 'a
20 val foldl
: ((elem
* 'a
) -> 'a
) -> 'a
-> vector
-> 'a
21 val foldr
: ((elem
* 'a
) -> 'a
) -> 'a
-> vector
-> 'a
25 Copyright (c
) Juan Jose Garcia Ripoll
.
28 Refer to the COPYRIGHT file for license conditions
33 Redistribution
and use
in source
and binary forms
, with or
34 without modification
, are permitted provided that the following
37 1. Redistributions
of source code must retain the above copyright
38 notice
, this list
of conditions
and the following disclaimer
.
40 2. Redistributions
in binary form must reproduce the above
41 copyright notice
, this list
of conditions
and the following
42 disclaimer
in the documentation
and/or other materials provided
43 with the distribution
.
45 3. All advertising materials mentioning features or use
of this
46 software must display the following acknowledgement
:
47 This product includes software developed by Juan Jose
50 4. The name
of Juan Jose Garcia Ripoll may not be used to endorse
51 or promote products derived from this software without
52 specific prior written permission
.
54 THIS SOFTWARE IS PROVIDED BY JUAN JOSE GARCIA RIPOLL ``AS IS
''
55 AND ANY EXPRESS OR IMPLIED WARRANTIES
, INCLUDING
, BUT NOT LIMITED
56 TO
, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
57 PARTICULAR PURPOSE ARE DISCLAIMED
. IN NO EVENT SHALL HE BE
58 LIABLE FOR ANY DIRECT
, INDIRECT
, INCIDENTAL
, SPECIAL
, EXEMPLARY
,
59 OR CONSEQUENTIAL
DAMAGES (INCLUDING
, BUT NOT LIMITED TO
,
60 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES
; LOSS OF USE
, DATA
,
61 OR PROFITS
; OR BUSINESS INTERRUPTION
) HOWEVER CAUSED AND ON ANY
62 THEORY OF LIABILITY
, WHETHER IN CONTRACT
, STRICT LIABILITY
, OR
63 TORT (INCLUDING NEGLIGENCE OR OTHERWISE
) ARISING IN ANY WAY OUT
64 OF THE USE OF THIS SOFTWARE
, EVEN IF ADVISED OF THE POSSIBILITY
71 val TIME
= ref (Time
.now())
74 (TIME
:= Time
.now(); ())
76 Time
.toMilliseconds(Time
.-(Time
.now(),!TIME
))
78 let val delta
= timerRead()
81 print (LargeInt
.toString delta
);
84 fun time f
= (timerOn(); f(); timerOff())
107 (f a
; app (a
+1, b
, f
))
111 fun app
' (a
, b
, d
, f
) =
113 (f a
; app
' (a
+d
, b
, d
, f
))
117 fun appi
' (a
, b
, d
, f
) =
119 (f a
; appi
' (a
+d
, b
, d
, f
))
126 Indices are a enumerable finite set
of data
with an order
and a map
127 to a continous nonnegative interval
of integers
. In the sample
128 implementation
, Index
, each index is a list
of integers
,
130 and each set
of indices is defined by a shape
, which has the same
131 shape
of an index but
with each integer incremented by one
135 type storage
= RowMajor | ColumnMajor
138 1) the underlying algorithms for this
structure
139 2) the most significant index
140 3) the index that varies more slowly
142 RowMajor means that first index is most significant
and varies
143 more slowly
, while ColumnMajor means that last index is the most
144 significant
and varies more slowly
. For instance
145 RowMajor
=> [0,0]<[0,1]<[1,0]<[1,1] (C
, C
++, Pascal
)
146 ColumnMajor
=> [0,0]>[1,0]>[0,1]>[1,1] (Fortran
)
149 Returns the last
/first index that belongs to the sed defined by
152 Checkes whether
'index
' belongs to the set defined by
'shape
'.
154 As we said
, indices can be sorted
and mapped to a finite set
of
155 integers
. 'toInt
' obtaines the integer number that corresponds to
158 It is equivalent to the partial evaluation
'toInt shape
' but
159 optimized for
'shape
'.
165 Obtain the following or previous index to the one we supply
.
166 next
and prev return an object
of type 'index option
' so that
167 if there is no such following
/previous
, the output is NONE
.
168 On the other hand
, next
'/prev
' raise an
exception when the
169 output is not well defined
and their output is always
of type
170 index
. next
/prev
/next
'/prev
' raise an
exception if 'index
'
171 does not belong to the set
of 'shape
'.
176 Iterates
'f
' over every index
of the set defined by
'shape
'.
177 'all
' stops when
'f
' first returns
false, 'any
' stops when
178 'f
' first returns
true and 'app
' does not stop
and discards the
182 Returns LESS
/GREATER
/EQUAL according to the total order which
183 is defined
in the set
of all indices
.
185 Reduced comparisons which are defined
in terms
of 'compare
'.
189 Checks whether
't
' conforms a valid shape or index
.
197 type indexer
= t
-> int
198 datatype storage
= RowMajor | ColumnMajor
204 val toInt
: t
-> t
-> int
205 val length
: t
-> int
208 val next
: t
-> t
-> t option
209 val prev
: t
-> t
-> t option
210 val next
' : t
-> t
-> t
211 val prev
' : t
-> t
-> t
212 val indexer
: t
-> (t
-> int)
214 val inBounds
: t
-> t
-> bool
215 val compare
: t
* t
-> order
216 val < : t
* t
-> bool
217 val > : t
* t
-> bool
218 val eq
: t
* t
-> bool
219 val <= : t
* t
-> bool
220 val >= : t
* t
-> bool
221 val <> : t
* t
-> bool
224 val validShape
: t
-> bool
225 val validIndex
: t
-> bool
227 val all
: t
-> (t
-> bool) -> bool
228 val any
: t
-> (t
-> bool) -> bool
229 val app
: t
-> (t
-> unit
) -> unit
231 structure Index
: INDEX
=
234 type indexer
= t
-> int
235 datatype storage
= RowMajor | ColumnMajor
240 val order
= ColumnMajor
242 fun validShape shape
= List.all (fn x
=> x
> 0) shape
244 fun validIndex index
= List.all (fn x
=> x
>= 0) index
246 fun toInt shape index
=
247 let fun loop ([], [], accum
, _
) = accum
248 |
loop ([], _
, _
, _
) = raise Index
249 |
loop (_
, [], _
, _
) = raise Index
250 |
loop (i
::ri
, l
::rl
, accum
, fac
) =
251 if (i
>= 0) andalso (i
< l
) then
252 loop (ri
, rl
, i
*fac
+ accum
, fac
*l
)
255 in loop (index
, shape
, 0, 1)
258 (* ----- CACHED LINEAR INDEXER
-----
260 An indexer is a function that takes a list
of
261 indices
, validates it
and produces a nonnegative
262 integer number
. In short
, the indexer is the
263 mapper from indices to element positions
in
266 'indexer
' builds such a mapper by optimizing
267 the most common cases
, which are
1d
and 2d
271 fun doindexer
[] _
= raise Shape
272 | doindexer
[a
] [dx
] =
273 let fun f
[x
] = if (x
> 0) andalso (x
< a
)
278 | doindexer
[a
,b
] [dx
, dy
] =
279 let fun f
[x
,y
] = if ((x
> 0) andalso (x
< a
) andalso
280 (y
> 0) andalso (y
< b
))
285 | doindexer
[a
,b
,c
] [dx
,dy
,dz
] =
286 let fun f
[x
,y
,z
] = if ((x
> 0) andalso (x
< a
) andalso
287 (y
> 0) andalso (y
< b
) andalso
288 (z
> 0) andalso (z
< c
))
289 then x
+ dy
* y
+ dz
* z
293 | doindexer shape memo
=
294 let fun f
[] [] accum
[] = accum
295 | f _ _ _
[] = raise Index
296 |
f (fac
::rf
) (ndx
::ri
) accum (dim
::rd
) =
297 if (ndx
>= 0) andalso (ndx
< dim
) then
298 f rf
ri (accum
+ ndx
* fac
) rd
305 let fun memoize accum
[] = []
306 | memoize
accum (dim
::rd
) =
307 accum
:: (memoize (dim
* accum
) rd
)
310 then doindexer
shape (memoize
1 shape
)
317 if b
< 0 then raise Shape
else a
* b
318 in foldl prod
1 shape
321 fun first shape
= map (fn x
=> 0) shape
324 |
last (size
:: rest
) =
327 else size
- 1 :: last rest
329 fun next
' [] [] = raise Subscript
330 | next
' _
[] = raise Index
331 | next
' [] _
= raise Index
332 | next
' (dimension
::restd
) (index
::resti
) =
333 if (index
+ 1) < dimension
334 then (index
+ 1) :: resti
335 else 0 :: (next
' restd resti
)
337 fun prev
' [] [] = raise Subscript
338 | prev
' _
[] = raise Index
339 | prev
' [] _
= raise Index
340 | prev
' (dimension
::restd
) (index
::resti
) =
342 then index
- 1 :: resti
343 else dimension
- 1 :: prev
' restd resti
345 fun next shape index
= (SOME (next
' shape index
)) handle
348 fun prev shape index
= (SOME (prev
' shape index
)) handle
351 fun inBounds shape index
=
352 ListPair.all (fn (x
,y
) => (x
>= 0) andalso (x
< y
))
355 fun compare ([],[]) = EQUAL
356 |
compare (_
, []) = raise Index
357 |
compare ([],_
) = raise Index
358 |
compare (a
::ra
, b
::rb
) =
359 case Int.compare (a
,b
) of
360 EQUAL
=> compare (ra
,rb
)
365 fun iterator a inner
=
366 let fun loop accum f
=
367 let fun innerloop i
=
369 then if inner (i
::accum
) f
377 fun build_iterator
[a
] =
378 let fun loop accum f
=
379 let fun innerloop i
=
389 |
build_iterator (a
::rest
) = iterator
a (build_iterator rest
)
391 fun all shape
= build_iterator shape
[]
395 fun iterator a inner
=
396 let fun loop accum f
=
397 let fun innerloop i
=
399 then if inner (i
::accum
) f
407 fun build_iterator
[a
] =
408 let fun loop accum f
=
409 let fun innerloop i
=
419 |
build_iterator (a
::rest
) = iterator
a (build_iterator rest
)
421 fun any shape
= build_iterator shape
[]
425 fun iterator a inner
=
426 let fun loop accum f
=
427 let fun innerloop i
=
429 then (inner (i
::accum
) f
;
436 fun build_iterator
[a
] =
437 let fun loop accum f
=
438 let fun innerloop i
=
440 then (f (i
::accum
); innerloop (i
+1))
446 |
build_iterator (a
::rest
) = iterator
a (build_iterator rest
)
448 fun app shape
= build_iterator shape
[]
451 fun a
< b
= compare(a
,b
) = LESS
452 fun a
> b
= compare(a
,b
) = GREATER
453 fun eq (a
, b
) = compare(a
,b
) = EQUAL
454 fun a
<> b
= not (a
= b
)
455 fun a
<= b
= not (a
> b
)
456 fun a
>= b
= not (a
< b
)
457 fun a
- b
= ListPair.map
Int.- (a
,b
)
461 Copyright (c
) Juan Jose Garcia Ripoll
.
464 Refer to the COPYRIGHT file for license conditions
470 Polymorphic tensors
of any
type. With
'tensor
' we denote
a (mutable
)
471 array
of any rank
, with as many indices
as one wishes
, and that may
472 be
traversed (map
, fold
, etc
) according to any
of those indices
.
475 Polymorphic tensor whose elements are all
of type 'a
.
476 val storage
= RowMajor | ColumnMajor
477 RowMajor
= data is stored
in consecutive cells
, first index
478 varying
fastest (FORTRAN convention
)
479 ColumnMajor
= data is stored
in consecutive cells
, last
480 index varying
fastest (C
,C
++,Pascal
,CommonLisp convention
)
481 new ([i1
,...,in],init
)
482 Build a new tensor
with n indices
, each
of sizes i1
...in,
484 fromArray (shape
,data
)
485 fromList (shape
,data
)
486 Use
'data
' to fill a tensor
of that shape
. An
exception is
487 raised
if 'data
' is too large or too small to properly
488 fill the vector
. Later use
of a
'data
' array is disregarded
489 -- one must think that the tensor now owns the array
.
493 Return the number
of elements
, the number
of indices
and
494 the
shape (size
of each index
) of the tensor
.
496 Return the data
of the tensor
in the form
of an array
.
497 Mutation
of this array may lead to unexpected behavior
.
499 sub (tensor
,[i1
,...,in])
500 update (tensor
,[i1
,...,in],new_value
)
501 Access the element that is indexed by the numbers
[i1
,..,in]
505 The same
as 'map
' and 'mapi
' but the function
'f
' outputs
506 nothing
and no new array is produced
, i
.e
. one only seeks
507 the side effect that
'f
' may produce
.
509 Apply function
'f
' to pairs
of elements
of 'a
' and 'b
'
510 and build a new tensor
with the output
. Both operands
511 must have the same shape or an
exception is raised
.
512 The procedure is sequential
, as specified by
'storage
'.
514 Fold
-left the elements
of tensor
'a
' along the n
-th
518 Folded boolean tests on the elements
of the tensor
.
523 structure Array
: ARRAY
524 structure Index
: INDEX
528 val new
: index
* 'a
-> 'a tensor
529 val tabulate
: index
* (index
-> 'a
) -> 'a tensor
530 val length
: 'a tensor
-> int
531 val rank
: 'a tensor
-> int
532 val shape
: 'a tensor
-> (index
)
533 val reshape
: index
-> 'a tensor
-> 'a tensor
534 val fromList
: index
* 'a list
-> 'a tensor
535 val fromArray
: index
* 'a array
-> 'a tensor
536 val toArray
: 'a tensor
-> 'a array
538 val sub
: 'a tensor
* index
-> 'a
539 val update
: 'a tensor
* index
* 'a
-> unit
540 val map
: ('a
-> 'b
) -> 'a tensor
-> 'b tensor
541 val map2
: ('a
* 'b
-> 'c
) -> 'a tensor
-> 'b tensor
-> 'c tensor
542 val app
: ('a
-> unit
) -> 'a tensor
-> unit
543 val appi
: (int * 'a
-> unit
) -> 'a tensor
-> unit
544 val foldl
: ('c
* 'a
-> 'c
) -> 'c
-> 'a tensor
-> int -> 'c tensor
545 val all
: ('a
-> bool) -> 'a tensor
-> bool
546 val any
: ('a
-> bool) -> 'a tensor
-> bool
550 Copyright (c
) Juan Jose Garcia Ripoll
.
553 Refer to the COPYRIGHT file for license conditions
556 structure Tensor
: TENSOR
=
558 structure Array
= Array
559 structure Index
= Index
562 type 'a tensor
= {shape
: index
, indexer
: Index
.indexer
, data
: 'a array
}
569 (*----- LOCALS
-----*)
571 fun make
' (shape
, data
) =
572 {shape
= shape
, indexer
= Index
.indexer shape
, data
= data
}
574 fun toInt
{shape
, indexer
, data
} index
= indexer index
577 let fun apply index
= f(Array
.sub(a
,index
)) in
578 Array
.tabulate(Array
.length a
, apply
)
581 fun splitList (l
as (a
::rest
), place
) =
582 let fun loop (left
,here
,right
) 0 = (List.rev left
,here
,right
)
583 |
loop (_
,_
,[]) place
= raise Index
584 |
loop (left
,here
,a
::right
) place
=
585 loop (here
::left
,a
,right
) (place
-1)
588 loop ([],a
,rest
) (List.length rest
- place
)
590 loop ([],a
,rest
) (place
- 1)
594 (*----- STRUCTURAL OPERATIONS
& QUERIES
------*)
596 fun new (shape
, init
) =
597 if not (Index
.validShape shape
) then
600 let val length
= Index
.length shape
in
602 indexer
= Index
.indexer shape
,
603 data
= Array
.array(length
,init
)}
606 fun toArray
{shape
, indexer
, data
} = data
608 fun length
{shape
, indexer
, data
} = Array
.length data
610 fun shape
{shape
, indexer
, data
} = shape
612 fun rank t
= List.length (shape t
)
614 fun reshape new_shape tensor
=
615 if Index
.validShape new_shape
then
616 case (Index
.length new_shape
) = length tensor
of
617 true => make
'(new_shape
, toArray tensor
)
618 |
false => raise Match
622 fun fromArray (s
, a
) =
623 case Index
.validShape s
andalso
624 ((Index
.length s
) = (Array
.length a
)) of
626 |
false => raise Shape
628 fun fromList (s
, a
) = fromArray (s
, Array
.fromList a
)
630 fun tabulate (shape
,f
) =
631 if Index
.validShape shape
then
632 let val last
= Index
.last shape
633 val length
= Index
.length shape
634 val c
= Array
.array(length
, f last
)
635 fun dotable (c
, indices
, i
) =
636 (Array
.update(c
, i
, f indices
);
639 | i
=> dotable(c
, Index
.prev
' shape indices
, i
-1))
641 make
'(shape
,dotable(c
, Index
.prev
' shape last
, length
-1))
646 (*----- ELEMENTWISE OPERATIONS
-----*)
648 fun sub (t
, index
) = Array
.sub(#data t
, toInt t index
)
650 fun update (t
, index
, value
) =
651 Array
.update(toArray t
, toInt t index
, value
)
653 fun map f
{shape
, indexer
, data
} =
654 {shape
= shape
, indexer
= indexer
, data
= array_map f data
}
657 let val {shape
, indexer
, data
} = t1
658 val {shape
=shape2
, indexer
=indexer2
, data
=data2
} = t2
659 fun apply i
= f (Array
.sub(data
,i
), Array
.sub(data2
,i
))
660 val len
= Array
.length data
662 if Index
.eq(shape
, shape2
) then
665 data
= Array
.tabulate(len
, apply
)}
670 fun appi f tensor
= Array
.appi
f (toArray tensor
)
672 fun app f tensor
= Array
.app
f (toArray tensor
)
675 let val a
= toArray tensor
676 in Loop
.all(0, length tensor
- 1, fn i
=>
681 let val a
= toArray tensor
682 in Loop
.any(0, length tensor
- 1, fn i
=>
686 fun foldl f init
{shape
, indexer
, data
=a
} index
=
687 let val (head
,lk
,tail
) = splitList(shape
, index
)
688 val li
= Index
.length head
689 val lj
= Index
.length tail
690 val c
= Array
.array(li
* lj
,init
)
691 fun loopi (0, _
, _
) = ()
692 |
loopi (i
, ia
, ic
) =
693 (Array
.update(c
, ic
, f(Array
.sub(c
,ic
), Array
.sub(a
,ia
)));
694 loopi (i
-1, ia
+1, ic
+1))
695 fun loopk (0, ia
, _
) = ia
696 |
loopk (k
, ia
, ic
) = (loopi (li
, ia
, ic
);
697 loopk (k
-1, ia
+li
, ic
))
698 fun loopj (0, _
, _
) = ()
699 |
loopj (j
, ia
, ic
) = loopj (j
-1, loopk(lk
,ia
,ic
), ic
+li
)
702 make
'(head @ tail
, c
)
709 Copyright (c
) Juan Jose Garcia Ripoll
.
712 Refer to the COPYRIGHT file for license conditions
716 MONO_TENSOR
- signature -
718 Monomorphic tensor
of arbitrary
data (not only numbers
). Operations
719 should be provided to run the data
in several ways
, according to one
723 The
type of the tensor itself
725 The
type of every element
726 val storage
= RowMajor | ColumnMajor
727 RowMajor
= data is stored
in consecutive cells
, first index
728 varying
fastest (FORTRAN convention
)
729 ColumnMajor
= data is stored
in consecutive cells
, last
730 index varying
fastest (C
,C
++,Pascal
,CommonLisp convention
)
731 new ([i1
,...,in],init
)
732 Build a new tensor
with n indices
, each
of sizes i1
...in,
734 fromArray (shape
,data
)
735 fromList (shape
,data
)
736 Use
'data
' to fill a tensor
of that shape
. An
exception is
737 raised
if 'data
' is too large or too small to properly
738 fill the vector
. Later use
of a
'data
' array is disregarded
739 -- one must think that the tensor now owns the array
.
743 Return the number
of elements
, the number
of indices
and
744 the
shape (size
of each index
) of the tensor
.
746 Return the data
of the tensor
in the form
of an array
.
747 Mutation
of this array may lead to unexpected behavior
.
748 The data
in the array is stored according to `storage
'.
750 sub (tensor
,[i1
,...,in])
751 update (tensor
,[i1
,...,in],new_value
)
752 Access the element that is indexed by the numbers
[i1
,..,in]
756 Produce a new array by mapping the function sequentially
757 as specified by
'storage
', to each element
of tensor
'a
'.
758 In
'mapi
' the function receives
a (indices
,value
) tuple
,
759 while in 'map
' it only receives the value
.
762 The same
as 'map
' and 'mapi
' but the function
'f
' outputs
763 nothing
and no new array is produced
, i
.e
. one only seeks
764 the side effect that
'f
' may produce
.
766 Apply function
'f
' to pairs
of elements
of 'a
' and 'b
'
767 and build a new tensor
with the output
. Both operands
768 must have the same shape or an
exception is raised
.
769 The procedure is sequential
, as specified by
'storage
'.
771 Fold
-left the elements
of tensor
'a
' along the n
-th
775 Folded boolean tests on the elements
of the tensor
.
778 Polymorphic versions
of map
, map2
, foldl
.
781 signature MONO_TENSOR
=
783 structure Array
: MONO_ARRAY
784 structure Index
: INDEX
790 val new
: index
* elem
-> tensor
791 val tabulate
: index
* (index
-> elem
) -> tensor
792 val length
: tensor
-> int
793 val rank
: tensor
-> int
794 val shape
: tensor
-> (index
)
795 val reshape
: index
-> tensor
-> tensor
796 val fromList
: index
* elem list
-> tensor
797 val fromArray
: index
* Array
.array
-> tensor
798 val toArray
: tensor
-> Array
.array
800 val sub
: tensor
* index
-> elem
801 val update
: tensor
* index
* elem
-> unit
802 val map
: (elem
-> elem
) -> tensor
-> tensor
803 val map2
: (elem
* elem
-> elem
) -> tensor
-> tensor
-> tensor
804 val app
: (elem
-> unit
) -> tensor
-> unit
805 val appi
: (int * elem
-> unit
) -> tensor
-> unit
806 val foldl
: (elem
* 'a
-> 'a
) -> 'a
-> tensor
-> tensor
807 val foldln
: (elem
* elem
-> elem
) -> elem
-> tensor
-> int -> tensor
808 val all
: (elem
-> bool) -> tensor
-> bool
809 val any
: (elem
-> bool) -> tensor
-> bool
811 val map
' : (elem
-> 'a
) -> tensor
-> 'a Tensor
.tensor
812 val map2
' : (elem
* elem
-> 'a
) -> tensor
-> tensor
-> 'a Tensor
.tensor
813 val foldl
' : ('a
* elem
-> 'a
) -> 'a
-> tensor
-> int -> 'a Tensor
.tensor
819 Guarantees a
structure with a minimal number
of mathematical operations
820 so
as to build an algebraic
structure named Tensor
.
833 val toString
: t
-> string
845 val *+ : t
* t
* t
-> t
846 val *- : t
* t
* t
-> t
847 val ** : t
* int -> t
853 val == : t
* t
-> bool
854 val != : t
* t
-> bool
856 val toString
: t
-> string
857 val fromInt
: int -> t
858 val scan
: (char
,'a
) StringCvt.reader
-> (t
,'a
) StringCvt.reader
861 signature INTEGRAL_NUMBER
=
865 val quot
: t
* t
-> t
870 val compare
: t
* t
-> order
871 val < : t
* t
-> bool
872 val > : t
* t
-> bool
873 val <= : t
* t
-> bool
874 val >= : t
* t
-> bool
880 signature FRACTIONAL_NUMBER
=
908 val atan2
: t
* t
-> t
911 signature REAL_NUMBER
=
913 include FRACTIONAL_NUMBER
915 val compare
: t
* t
-> order
916 val < : t
* t
-> bool
917 val > : t
* t
-> bool
918 val <= : t
* t
-> bool
919 val >= : t
* t
-> bool
925 signature COMPLEX_NUMBER
=
927 include FRACTIONAL_NUMBER
929 structure Real : REAL_NUMBER
932 val make
: real * real -> t
933 val split
: t
-> real * real
934 val realPart
: t
-> real
935 val imagPart
: t
-> real
939 structure INumber
: INTEGRAL_NUMBER
=
951 let val x
= loop (Int.div(n
, 2))
952 val m
= Int.mod(n
, 2)
964 fun signum i
= case compare(i
, 0) of
972 fun a
!= b
= (a
<> b
)
973 fun *+(b
,c
,a
) = b
* c
+ a
974 fun *-(b
,c
,a
) = b
* c
- b
976 fun scan getc
= Int.scan
StringCvt.DEC getc
979 structure RNumber
: REAL_NUMBER
=
987 fun signum x
= case compare(x
,0.0) of
992 fun recip x
= 1.0 / x
999 let val x
= loop (Int.div(n
, 2))
1000 val m
= Int.mod(n
, 2)
1012 fun max (a
, b
) = if a
< b
then b
else a
1013 fun min (a
, b
) = if a
< b
then a
else b
1015 fun asinh x
= ln (x
+ sqrt(1.0 + x
* x
))
1016 fun acosh x
= ln (x
+ (x
+ 1.0) * sqrt((x
- 1.0)/(x
+ 1.0)))
1017 fun atanh x
= ln ((1.0 + x
) / sqrt(1.0 - x
* x
))
1021 Complex(R
) - Functor
-
1023 Provides support for complex numbers based on tuples
. Should be
1024 highly efficient
as most operations can be inlined
.
1027 structure CNumber
: COMPLEX_NUMBER
=
1029 structure Real = RNumber
1031 type t
= Real.t
* Real.t
1034 val zero
= (0.0,0.0)
1036 val pi
= (Real.pi
, 0.0)
1037 val e
= (Real.e
, 0.0)
1039 fun make (r
,i
) = (r
,i
) : t
1041 fun realPart (r
,_
) = r
1042 fun imagPart (_
,i
) = i
1044 fun abs2 (r
,i
) = Real.+(Real.*(r
,r
),Real.*(i
,i
)) (* FIXME
!!! *)
1045 fun arg (r
,i
) = Real.atan2(i
,r
)
1046 fun modulus z
= Real.sqrt(abs2 z
)
1047 fun abs z
= (modulus z
, 0.0)
1048 fun signum (z
as (r
,i
)) =
1049 let val m
= modulus z
1050 in (Real./(r
,m
), Real./(i
,m
))
1053 fun ~
(r1
,i1
) = (Real.~ r1
, Real.~ i1
)
1054 fun (r1
,i1
) + (r2
,i2
) = (Real.+(r1
,r2
), Real.+(i1
,i2
))
1055 fun (r1
,i1
) - (r2
,i2
) = (Real.-(r1
,r2
), Real.-(i1
,i1
))
1056 fun (r1
,i1
) * (r2
,i2
) = (Real.-(Real.*(r1
,r2
),Real.*(i1
,i2
)),
1057 Real.+(Real.*(r1
,i2
),Real.*(r2
,i1
)))
1058 fun (r1
,i1
) / (r2
,i2
) =
1059 let val modulus
= abs2(r2
,i2
)
1060 val (nr
,ni
) = (r1
,i1
) * (r2
,i2
)
1062 (Real./(nr
,modulus
), Real./(ni
,modulus
))
1064 fun *+((r1
,i1
),(r2
,i2
),(r0
,i0
)) =
1065 (Real.*+(Real.~ i1
, i2
, Real.*+(r1
,r2
,r0
)),
1066 Real.*+(r2
, i2
, Real.*+(r1
,i2
,i0
)))
1067 fun *-((r1
,i1
),(r2
,i2
),(r0
,i0
)) =
1068 (Real.*+(Real.~ i1
, i2
, Real.*-(r1
,r2
,r0
)),
1069 Real.*+(r2
, i2
, Real.*-(r1
,i2
,i0
)))
1073 let fun loop
0 = one
1076 let val x
= loop (Int.div(n
, 2))
1077 val m
= Int.mod(n
, 2)
1089 fun recip (r1
, i1
) =
1090 let val modulus
= abs2(r1
, i1
)
1091 in (Real./(r1
, modulus
), Real./(Real.~ i1
, modulus
))
1093 fun ==(z
, w
) = Real.==(realPart z
, realPart w
) andalso Real.==(imagPart z
, imagPart w
)
1094 fun !=(z
, w
) = Real.!=(realPart z
, realPart w
) andalso Real.!=(imagPart z
, imagPart w
)
1095 fun fromInt i
= (Real.fromInt i
, 0.0)
1096 fun toString (r
,i
) =
1097 String.concat
["(",Real.toString r
,",",Real.toString i
,")"]
1100 let val expx
= Real.exp x
1101 in (Real.*(x
, (Real.cos y
)), Real.*(x
, (Real.sin y
)))
1105 val half
= Real.recip (Real.fromInt
2)
1107 fun sqrt (z
as (x
,y
)) =
1108 if Real.==(x
, 0.0) andalso Real.==(y
, 0.0) then
1111 let val m
= Real.+(modulus z
, Real.abs x
)
1112 val u
' = Real.sqrt (Real.*(m
, half
))
1113 val v
' = Real./(Real.abs y
, Real.+(u
',u
'))
1114 val (u
,v
) = if Real.<(x
, 0.0) then (v
',u
') else (u
',v
')
1115 in (u
, if Real.<(y
, 0.0) then Real.~ v
else v
)
1118 fun ln z
= (Real.ln (modulus z
), arg z
)
1125 fun sin (x
, y
) = (Real.*(Real.sin x
, Real.cosh y
),
1126 Real.*(Real.cos x
, Real.sinh y
))
1127 fun cos (x
, y
) = (Real.*(Real.cos x
, Real.cosh y
),
1128 Real.~
(Real.*(Real.sin x
, Real.sinh y
)))
1130 let val (sx
, cx
) = (Real.sin x
, Real.cos x
)
1131 val (shy
, chy
) = (Real.sinh y
, Real.cosh y
)
1132 val a
= (Real.*(sx
, chy
), Real.*(cx
, shy
))
1133 val b
= (Real.*(cx
, chy
), Real.*(Real.~ sx
, shy
))
1137 fun sinh (x
, y
) = (Real.*(Real.cos y
, Real.sinh x
),
1138 Real.*(Real.sin y
, Real.cosh x
))
1139 fun cosh (x
, y
) = (Real.*(Real.cos y
, Real.cosh x
),
1140 Real.*(Real.sin y
, Real.sinh x
))
1142 let val (sy
, cy
) = (Real.sin y
, Real.cos y
)
1143 val (shx
, chx
) = (Real.sinh x
, Real.cosh x
)
1144 val a
= (Real.*(cy
, shx
), Real.*(sy
, chx
))
1145 val b
= (Real.*(cy
, chx
), Real.*(sy
, shx
))
1149 fun asin (z
as (x
,y
)) =
1150 let val w
= sqrt (one
- z
* z
)
1151 val (x
',y
') = ln ((Real.~ y
, x
) + w
)
1155 fun acos (z
as (x
,y
)) =
1156 let val (x
', y
') = sqrt (one
+ z
* z
)
1157 val (x
'', y
'') = ln (z
+ (Real.~ y
', x
'))
1158 in (y
'', Real.~ x
'')
1161 fun atan (z
as (x
,y
)) =
1162 let val w
= sqrt (one
+ z
*z
)
1163 val (x
',y
') = ln ((Real.-(1.0, y
), x
) / w
)
1167 fun atan2 (y
, x
) = atan(y
/ x
)
1169 fun asinh x
= ln (x
+ sqrt(one
+ x
* x
))
1170 fun acosh x
= ln (x
+ (x
+ one
) * sqrt((x
- one
)/(x
+ one
)))
1171 fun atanh x
= ln ((one
+ x
) / sqrt(one
- x
* x
))
1174 let val scanner
= Real.scan getc
1176 case scanner stream
of
1179 case scanner rest
of
1181 |
SOME (b
, rest
) => SOME (make(a
,b
), rest
)
1184 end (* ComplexNumber
*)
1187 Copyright (c
) Juan Jose Garcia Ripoll
.
1188 All rights reserved
.
1190 Refer to the COPYRIGHT file for license conditions
1193 structure INumberArray
=
1196 type array
= INumber
.t array
1197 type vector
= INumber
.t vector
1198 type elem
= INumber
.t
1202 type vector
= INumber
.t
Vector.vector
1203 type elem
= INumber
.t
1205 fun map f a
= tabulate(length a
, fn x
=> (f (sub(a
,x
))))
1206 fun mapi f a
= tabulate(length a
, fn x
=> (f (x
,sub(a
,x
))))
1207 fun map2 f a b
= tabulate(length a
, fn x
=> (f(sub(a
,x
),sub(b
,x
))))
1210 structure RNumberArray
=
1213 val sub
= Unsafe
.Real64Array
.sub
1214 val update
= Unsafe
.Real64Array
.update
1215 fun map f a
= tabulate(length a
, fn x
=> (f (sub(a
,x
))))
1216 fun mapi f a
= tabulate(length a
, fn x
=> (f (x
,sub(a
,x
))))
1217 fun map2 f a b
= tabulate(length a
, fn x
=> (f(sub(a
,x
),sub(b
,x
))))
1220 (*--------------------- COMPLEX ARRAY
-------------------------*)
1222 structure BasicCNumberArray
=
1224 structure Complex
: COMPLEX_NUMBER
= CNumber
1225 structure Array
: MONO_ARRAY
= RNumberArray
1227 type elem
= Complex
.t
1228 type array
= Array
.array
* Array
.array
1230 val maxLen
= Array
.maxLen
1232 fun length (a
,b
) = Array
.length a
1234 fun sub ((a
,b
),index
) = Complex
.make(Array
.sub(a
,index
),Array
.sub(b
,index
))
1236 fun update ((a
,b
),index
,z
) =
1237 let val (re
,im
) = Complex
.split z
in
1238 Array
.update(a
, index
, re
);
1239 Array
.update(b
, index
, im
)
1243 fun makeRange (a
, start
, NONE
) = makeRange(a
, start
, SOME (length a
- 1))
1244 |
makeRange (a
, start
, SOME last
) =
1245 let val len
= length a
1246 val diff
= last
- start
1248 if (start
>= len
) orelse (last
>= len
) then
1250 else if diff
< 0 then
1253 (a
, start
, diff
+ 1)
1258 fun array (size
,z
:elem
) =
1259 let val realsize
= size
* 2
1260 val r
= Complex
.realPart z
1261 val i
= Complex
.imagPart z
in
1262 (Array
.array(size
,r
), Array
.array(size
,i
))
1265 fun zeroarray size
=
1266 (Array
.array(size
,Complex
.Real.zero
),
1267 Array
.array(size
,Complex
.Real.zero
))
1269 fun tabulate (size
,f
) =
1270 let val a
= array(size
, Complex
.zero
)
1274 |
false => (update(a
, i
, f i
); loop (i
+1))
1280 let val length
= List.length list
1281 val a
= zeroarray length
1282 fun loop (_
, []) = a
1283 |
loop (i
, z
::rest
) = (update(a
, i
, z
);
1290 let val (a
, start
, len
) = makeRange range
1291 fun copy i
= sub(a
, i
+ start
)
1292 in tabulate(len
, copy
)
1295 fun concat array_list
=
1296 let val total_length
= foldl (op +) 0 (map length array_list
)
1297 val a
= array(total_length
, Complex
.zero
)
1298 fun copy (_
, []) = a
1299 |
copy (pos
, v
::rest
) =
1303 |
false => (update(a
, i
+pos
, sub(v
, i
)); loop (i
-1))
1304 in (loop (length v
- 1); copy(length v
+ pos
, rest
))
1310 fun copy
{src
: array
, si
: int, len
: int option
, dst
: array
, di
: int } =
1311 let val (a
, ia
, la
) = makeRange (src
, si
, len
)
1312 val (b
, ib
, lb
) = makeRange (dst
, di
, len
)
1316 |
false => (update(b
, i
+ib
, sub(a
, i
+ia
)); copy (i
-1))
1322 fun modifyi f range
=
1323 let val (a
, start
, len
) = makeRange range
1324 val last
= start
+ len
1328 |
false => (update(a
, i
, f(i
, sub(a
,i
))); loop (i
+1))
1333 let val last
= length a
1337 |
false => (update(a
, i
, f(sub(a
,i
))); loop (i
+1))
1342 let val size
= length a
1346 |
false => (f(sub(a
,i
)); loop (i
+1))
1352 let val (a
, start
, len
) = makeRange range
1353 val last
= start
+ len
1357 |
false => (f(i
, sub(a
,i
)); loop (i
+1))
1363 let val len
= length a
1364 val c
= zeroarray len
1366 | loop i
= (update(a
, i
, f(sub(a
,i
))); loop (i
-1))
1371 let val len
= length a
1372 val c
= zeroarray len
1374 | loop i
= (update(c
, i
, f(sub(a
,i
),sub(b
,i
)));
1380 let val (a
, start
, len
) = makeRange range
1381 fun rule i
= f (i
+start
, sub(a
, i
+start
))
1382 in tabulate(len
, rule
)
1385 fun foldli f init range
=
1386 let val (a
, start
, len
) = makeRange range
1387 val last
= start
+ len
- 1
1388 fun loop (i
, accum
) =
1391 |
false => loop (i
+1, f(i
, sub(a
,i
), accum
))
1392 in loop (start
, init
)
1395 fun foldri f init range
=
1396 let val (a
, start
, len
) = makeRange range
1397 val last
= start
+ len
- 1
1398 fun loop (i
, accum
) =
1401 |
false => loop (i
-1, f(i
, sub(a
,i
), accum
))
1402 in loop (last
, init
)
1405 fun foldl f init a
= foldli (fn (_
, a
, x
) => f(a
,x
)) init (a
,0,NONE
)
1406 fun foldr f init a
= foldri (fn (_
, x
, a
) => f(x
,a
)) init (a
,0,NONE
)
1408 end (* BasicCNumberArray
*)
1411 structure CNumberArray
=
1415 open BasicCNumberArray
1418 type vector
= Vector.vector
1419 open BasicCNumberArray
1420 end (* CNumberArray
*)
1421 structure INumber
: INTEGRAL_NUMBER
=
1432 let val x
= loop (Int.div(n
, 2))
1433 val m
= Int.mod(n
, 2)
1444 fun signum i
= case compare(i
, 0) of
1451 fun a
!= b
= (a
<> b
)
1452 fun *+(b
,c
,a
) = b
* c
+ a
1453 fun *-(b
,c
,a
) = b
* c
- b
1454 fun scan getc
= Int.scan
StringCvt.DEC getc
1456 structure RNumber
: REAL_NUMBER
=
1463 fun signum x
= case compare(x
,0.0) of
1467 fun recip x
= 1.0 / x
1470 let fun loop
0 = one
1473 let val x
= loop (Int.div(n
, 2))
1474 val m
= Int.mod(n
, 2)
1485 fun max (a
, b
) = if a
< b
then b
else a
1486 fun min (a
, b
) = if a
< b
then a
else b
1487 fun asinh x
= ln (x
+ sqrt(1.0 + x
* x
))
1488 fun acosh x
= ln (x
+ (x
+ 1.0) * sqrt((x
- 1.0)/(x
+ 1.0)))
1489 fun atanh x
= ln ((1.0 + x
) / sqrt(1.0 - x
* x
))
1492 Complex(R
) - Functor
-
1493 Provides support for complex numbers based on tuples
. Should be
1494 highly efficient
as most operations can be inlined
.
1496 structure CNumber
: COMPLEX_NUMBER
=
1498 structure Real = RNumber
1499 type t
= Real.t
* Real.t
1501 val zero
= (0.0,0.0)
1503 val pi
= (Real.pi
, 0.0)
1504 val e
= (Real.e
, 0.0)
1505 fun make (r
,i
) = (r
,i
) : t
1507 fun realPart (r
,_
) = r
1508 fun imagPart (_
,i
) = i
1509 fun abs2 (r
,i
) = Real.+(Real.*(r
,r
),Real.*(i
,i
)) (* FIXME
!!! *)
1510 fun arg (r
,i
) = Real.atan2(i
,r
)
1511 fun modulus z
= Real.sqrt(abs2 z
)
1512 fun abs z
= (modulus z
, 0.0)
1513 fun signum (z
as (r
,i
)) =
1514 let val m
= modulus z
1515 in (Real./(r
,m
), Real./(i
,m
))
1517 fun ~
(r1
,i1
) = (Real.~ r1
, Real.~ i1
)
1518 fun (r1
,i1
) + (r2
,i2
) = (Real.+(r1
,r2
), Real.+(i1
,i2
))
1519 fun (r1
,i1
) - (r2
,i2
) = (Real.-(r1
,r2
), Real.-(i1
,i1
))
1520 fun (r1
,i1
) * (r2
,i2
) = (Real.-(Real.*(r1
,r2
),Real.*(i1
,i2
)),
1521 Real.+(Real.*(r1
,i2
),Real.*(r2
,i1
)))
1522 fun (r1
,i1
) / (r2
,i2
) =
1523 let val modulus
= abs2(r2
,i2
)
1524 val (nr
,ni
) = (r1
,i1
) * (r2
,i2
)
1526 (Real./(nr
,modulus
), Real./(ni
,modulus
))
1528 fun *+((r1
,i1
),(r2
,i2
),(r0
,i0
)) =
1529 (Real.*+(Real.~ i1
, i2
, Real.*+(r1
,r2
,r0
)),
1530 Real.*+(r2
, i2
, Real.*+(r1
,i2
,i0
)))
1531 fun *-((r1
,i1
),(r2
,i2
),(r0
,i0
)) =
1532 (Real.*+(Real.~ i1
, i2
, Real.*-(r1
,r2
,r0
)),
1533 Real.*+(r2
, i2
, Real.*-(r1
,i2
,i0
)))
1536 let fun loop
0 = one
1539 let val x
= loop (Int.div(n
, 2))
1540 val m
= Int.mod(n
, 2)
1551 fun recip (r1
, i1
) =
1552 let val modulus
= abs2(r1
, i1
)
1553 in (Real./(r1
, modulus
), Real./(Real.~ i1
, modulus
))
1555 fun ==(z
, w
) = Real.==(realPart z
, realPart w
) andalso Real.==(imagPart z
, imagPart w
)
1556 fun !=(z
, w
) = Real.!=(realPart z
, realPart w
) andalso Real.!=(imagPart z
, imagPart w
)
1557 fun fromInt i
= (Real.fromInt i
, 0.0)
1558 fun toString (r
,i
) =
1559 String.concat
["(",Real.toString r
,",",Real.toString i
,")"]
1561 let val expx
= Real.exp x
1562 in (Real.*(x
, (Real.cos y
)), Real.*(x
, (Real.sin y
)))
1565 val half
= Real.recip (Real.fromInt
2)
1567 fun sqrt (z
as (x
,y
)) =
1568 if Real.==(x
, 0.0) andalso Real.==(y
, 0.0) then
1571 let val m
= Real.+(modulus z
, Real.abs x
)
1572 val u
' = Real.sqrt (Real.*(m
, half
))
1573 val v
' = Real./(Real.abs y
, Real.+(u
',u
'))
1574 val (u
,v
) = if Real.<(x
, 0.0) then (v
',u
') else (u
',v
')
1575 in (u
, if Real.<(y
, 0.0) then Real.~ v
else v
)
1578 fun ln z
= (Real.ln (modulus z
), arg z
)
1583 fun sin (x
, y
) = (Real.*(Real.sin x
, Real.cosh y
),
1584 Real.*(Real.cos x
, Real.sinh y
))
1585 fun cos (x
, y
) = (Real.*(Real.cos x
, Real.cosh y
),
1586 Real.~
(Real.*(Real.sin x
, Real.sinh y
)))
1588 let val (sx
, cx
) = (Real.sin x
, Real.cos x
)
1589 val (shy
, chy
) = (Real.sinh y
, Real.cosh y
)
1590 val a
= (Real.*(sx
, chy
), Real.*(cx
, shy
))
1591 val b
= (Real.*(cx
, chy
), Real.*(Real.~ sx
, shy
))
1594 fun sinh (x
, y
) = (Real.*(Real.cos y
, Real.sinh x
),
1595 Real.*(Real.sin y
, Real.cosh x
))
1596 fun cosh (x
, y
) = (Real.*(Real.cos y
, Real.cosh x
),
1597 Real.*(Real.sin y
, Real.sinh x
))
1599 let val (sy
, cy
) = (Real.sin y
, Real.cos y
)
1600 val (shx
, chx
) = (Real.sinh x
, Real.cosh x
)
1601 val a
= (Real.*(cy
, shx
), Real.*(sy
, chx
))
1602 val b
= (Real.*(cy
, chx
), Real.*(sy
, shx
))
1605 fun asin (z
as (x
,y
)) =
1606 let val w
= sqrt (one
- z
* z
)
1607 val (x
',y
') = ln ((Real.~ y
, x
) + w
)
1610 fun acos (z
as (x
,y
)) =
1611 let val (x
', y
') = sqrt (one
+ z
* z
)
1612 val (x
'', y
'') = ln (z
+ (Real.~ y
', x
'))
1613 in (y
'', Real.~ x
'')
1615 fun atan (z
as (x
,y
)) =
1616 let val w
= sqrt (one
+ z
*z
)
1617 val (x
',y
') = ln ((Real.-(1.0, y
), x
) / w
)
1620 fun atan2 (y
, x
) = atan(y
/ x
)
1621 fun asinh x
= ln (x
+ sqrt(one
+ x
* x
))
1622 fun acosh x
= ln (x
+ (x
+ one
) * sqrt((x
- one
)/(x
+ one
)))
1623 fun atanh x
= ln ((one
+ x
) / sqrt(one
- x
* x
))
1625 let val scanner
= Real.scan getc
1627 case scanner stream
of
1630 case scanner rest
of
1632 |
SOME (b
, rest
) => SOME (make(a
,b
), rest
)
1634 end (* ComplexNumber
*)
1636 Copyright (c
) Juan Jose Garcia Ripoll
.
1637 All rights reserved
.
1638 Refer to the COPYRIGHT file for license conditions
1640 structure PrettyPrint
:>
1645 Complex
of CNumber
.t |
1647 val list
: ('a
-> string) -> 'a list
-> unit
1648 val intList
: int list
-> unit
1649 val realList
: real list
-> unit
1650 val stringList
: string list
-> unit
1651 val array
: ('a
-> string) -> 'a array
-> unit
1652 val intArray
: int array
-> unit
1653 val realArray
: real array
-> unit
1654 val stringArray
: string array
-> unit
1656 int -> ((int * 'a
-> unit
) -> 'b
-> unit
) -> ('a
-> string) -> 'b
-> unit
1657 val print
: modifier list
-> unit
1663 Complex
of CNumber
.t |
1665 fun list _
[] = print
"[]"
1666 | list
cvt (a
::resta
) =
1667 let fun loop a
[] = (print(cvt a
); print
"]")
1668 | loop
a (b
::restb
) = (print(cvt a
); print
", "; loop b restb
)
1673 fun boolList a
= list
Bool.toString a
1674 fun intList a
= list
Int.toString a
1675 fun realList a
= list
Real.toString a
1676 fun stringList a
= list (fn x
=> x
) a
1678 let val length
= Array
.length a
- 1
1679 fun print_one (i
,x
) =
1680 (print(cvt x
); if not(i
= length
) then print
", " else ())
1682 Array
.appi print_one a
1684 fun boolArray a
= array
Bool.toString a
1685 fun intArray a
= array
Int.toString a
1686 fun realArray a
= array
Real.toString a
1687 fun stringArray a
= array (fn x
=> x
) a
1688 fun sequence length appi cvt seq
=
1689 let val length
= length
- 1
1690 fun print_one (i
:int,x
) =
1691 (print(cvt x
); if not(i
= length
) then print
", " else ())
1698 let fun printer (Int a
) = INumber
.toString a
1699 |
printer (Real a
) = RNumber
.toString a
1700 |
printer (Complex a
) = CNumber
.toString a
1701 |
printer (String a
) = a
1702 in List.app (fn x
=> (TextIO.print (printer x
))) b
1704 end (* PrettyPrint
*)
1705 fun print
' x
= List.app print x
1707 Copyright (c
) Juan Jose Garcia Ripoll
.
1708 All rights reserved
.
1709 Refer to the COPYRIGHT file for license conditions
1711 structure INumberArray
=
1714 type array
= INumber
.t array
1715 type vector
= INumber
.t vector
1716 type elem
= INumber
.t
1720 type vector
= INumber
.t
Vector.vector
1721 type elem
= INumber
.t
1723 fun map f a
= tabulate(length a
, fn x
=> (f (sub(a
,x
))))
1724 fun mapi f a
= tabulate(length a
, fn x
=> (f (x
,sub(a
,x
))))
1725 fun map2 f a b
= tabulate(length a
, fn x
=> (f(sub(a
,x
),sub(b
,x
))))
1727 structure RNumberArray
=
1730 val sub
= Unsafe
.Real64Array
.sub
1731 val update
= Unsafe
.Real64Array
.update
1732 fun map f a
= tabulate(length a
, fn x
=> (f (sub(a
,x
))))
1733 fun mapi f a
= tabulate(length a
, fn x
=> (f (x
,sub(a
,x
))))
1734 fun map2 f a b
= tabulate(length a
, fn x
=> (f(sub(a
,x
),sub(b
,x
))))
1736 (*--------------------- COMPLEX ARRAY
-------------------------*)
1737 structure BasicCNumberArray
=
1739 structure Complex
: COMPLEX_NUMBER
= CNumber
1740 structure Array
: MONO_ARRAY
= RNumberArray
1741 type elem
= Complex
.t
1742 type array
= Array
.array
* Array
.array
1743 val maxLen
= Array
.maxLen
1744 fun length (a
,b
) = Array
.length a
1745 fun sub ((a
,b
),index
) = Complex
.make(Array
.sub(a
,index
),Array
.sub(b
,index
))
1746 fun update ((a
,b
),index
,z
) =
1747 let val (re
,im
) = Complex
.split z
in
1748 Array
.update(a
, index
, re
);
1749 Array
.update(b
, index
, im
)
1752 fun makeRange (a
, start
, NONE
) = makeRange(a
, start
, SOME (length a
- 1))
1753 |
makeRange (a
, start
, SOME last
) =
1754 let val len
= length a
1755 val diff
= last
- start
1757 if (start
>= len
) orelse (last
>= len
) then
1759 else if diff
< 0 then
1762 (a
, start
, diff
+ 1)
1765 fun array (size
,z
:elem
) =
1766 let val realsize
= size
* 2
1767 val r
= Complex
.realPart z
1768 val i
= Complex
.imagPart z
in
1769 (Array
.array(size
,r
), Array
.array(size
,i
))
1771 fun zeroarray size
=
1772 (Array
.array(size
,Complex
.Real.zero
),
1773 Array
.array(size
,Complex
.Real.zero
))
1774 fun tabulate (size
,f
) =
1775 let val a
= array(size
, Complex
.zero
)
1779 |
false => (update(a
, i
, f i
); loop (i
+1))
1784 let val length
= List.length list
1785 val a
= zeroarray length
1786 fun loop (_
, []) = a
1787 |
loop (i
, z
::rest
) = (update(a
, i
, z
);
1793 let val (a
, start
, len
) = makeRange range
1794 fun copy i
= sub(a
, i
+ start
)
1795 in tabulate(len
, copy
)
1797 fun concat array_list
=
1798 let val total_length
= foldl (op +) 0 (map length array_list
)
1799 val a
= array(total_length
, Complex
.zero
)
1800 fun copy (_
, []) = a
1801 |
copy (pos
, v
::rest
) =
1805 |
false => (update(a
, i
+pos
, sub(v
, i
)); loop (i
-1))
1806 in (loop (length v
- 1); copy(length v
+ pos
, rest
))
1811 fun copy
{src
: array
, si
: int, len
: int option
, dst
: array
, di
: int } =
1812 let val (a
, ia
, la
) = makeRange (src
, si
, len
)
1813 val (b
, ib
, lb
) = makeRange (dst
, di
, len
)
1817 |
false => (update(b
, i
+ib
, sub(a
, i
+ia
)); copy (i
-1))
1821 fun modifyi f range
=
1822 let val (a
, start
, len
) = makeRange range
1823 val last
= start
+ len
1827 |
false => (update(a
, i
, f(i
, sub(a
,i
))); loop (i
+1))
1831 let val last
= length a
1835 |
false => (update(a
, i
, f(sub(a
,i
))); loop (i
+1))
1839 let val size
= length a
1843 |
false => (f(sub(a
,i
)); loop (i
+1))
1848 let val (a
, start
, len
) = makeRange range
1849 val last
= start
+ len
1853 |
false => (f(i
, sub(a
,i
)); loop (i
+1))
1858 let val len
= length a
1859 val c
= zeroarray len
1861 | loop i
= (update(a
, i
, f(sub(a
,i
))); loop (i
-1))
1865 let val len
= length a
1866 val c
= zeroarray len
1868 | loop i
= (update(c
, i
, f(sub(a
,i
),sub(b
,i
)));
1873 let val (a
, start
, len
) = makeRange range
1874 fun rule i
= f (i
+start
, sub(a
, i
+start
))
1875 in tabulate(len
, rule
)
1877 fun foldli f init range
=
1878 let val (a
, start
, len
) = makeRange range
1879 val last
= start
+ len
- 1
1880 fun loop (i
, accum
) =
1883 |
false => loop (i
+1, f(i
, sub(a
,i
), accum
))
1884 in loop (start
, init
)
1886 fun foldri f init range
=
1887 let val (a
, start
, len
) = makeRange range
1888 val last
= start
+ len
- 1
1889 fun loop (i
, accum
) =
1892 |
false => loop (i
-1, f(i
, sub(a
,i
), accum
))
1893 in loop (last
, init
)
1895 fun foldl f init a
= foldli (fn (_
, a
, x
) => f(a
,x
)) init (a
,0,NONE
)
1896 fun foldr f init a
= foldri (fn (_
, x
, a
) => f(x
,a
)) init (a
,0,NONE
)
1898 end (* BasicCNumberArray
*)
1899 structure CNumberArray
=
1903 open BasicCNumberArray
1906 type vector
= Vector.vector
1907 open BasicCNumberArray
1908 end (* CNumberArray
*)
1911 structure Number
= INumber
1912 structure Array
= INumberArray
1914 Copyright (c
) Juan Jose Garcia Ripoll
.
1915 All rights reserved
.
1916 Refer to the COPYRIGHT file for license conditions
1918 structure MonoTensor
=
1921 structure Array
= Array
1923 structure Index
= Index
1924 type elem
= Array
.elem
1925 type index
= Index
.t
1926 type tensor
= {shape
: index
, indexer
: Index
.indexer
, data
: Array
.array
}
1932 (*----- LOCALS
-----*)
1933 fun make
' (shape
, data
) =
1934 {shape
= shape
, indexer
= Index
.indexer shape
, data
= data
}
1935 fun toInt
{shape
, indexer
, data
} index
= indexer index
1936 fun splitList (l
as (a
::rest
), place
) =
1937 let fun loop (left
,here
,right
) 0 = (List.rev left
,here
,right
)
1938 |
loop (_
,_
,[]) place
= raise Index
1939 |
loop (left
,here
,a
::right
) place
=
1940 loop (here
::left
,a
,right
) (place
-1)
1943 loop ([],a
,rest
) (List.length rest
- place
)
1945 loop ([],a
,rest
) (place
- 1)
1948 (*----- STRUCTURAL OPERATIONS
& QUERIES
------*)
1949 fun new (shape
, init
) =
1950 if not (Index
.validShape shape
) then
1953 let val length
= Index
.length shape
in
1955 indexer
= Index
.indexer shape
,
1956 data
= Array
.array(length
,init
)}
1958 fun toArray
{shape
, indexer
, data
} = data
1959 fun length
{shape
, indexer
, data
} = Array
.length data
1960 fun shape
{shape
, indexer
, data
} = shape
1961 fun rank t
= List.length (shape t
)
1962 fun reshape new_shape tensor
=
1963 if Index
.validShape new_shape
then
1964 case (Index
.length new_shape
) = length tensor
of
1965 true => make
'(new_shape
, toArray tensor
)
1966 |
false => raise Match
1969 fun fromArray (s
, a
) =
1970 case Index
.validShape s
andalso
1971 ((Index
.length s
) = (Array
.length a
)) of
1973 |
false => raise Shape
1974 fun fromList (s
, a
) = fromArray (s
, Array
.fromList a
)
1975 fun tabulate (shape
,f
) =
1976 if Index
.validShape shape
then
1977 let val last
= Index
.last shape
1978 val length
= Index
.length shape
1979 val c
= Array
.array(length
, f last
)
1980 fun dotable (c
, indices
, i
) =
1981 (Array
.update(c
, i
, f indices
);
1984 else dotable(c
, Index
.prev
' shape indices
, i
-1))
1985 in make
'(shape
,dotable(c
, Index
.prev
' shape last
, length
-2))
1989 (*----- ELEMENTWISE OPERATIONS
-----*)
1990 fun sub (t
, index
) = Array
.sub(#data t
, toInt t index
)
1991 fun update (t
, index
, value
) =
1992 Array
.update(toArray t
, toInt t index
, value
)
1993 fun map f
{shape
, indexer
, data
} =
1994 {shape
= shape
, indexer
= indexer
, data
= Array
.map f data
}
1996 let val {shape
=shape1
, indexer
=indexer1
, data
=data1
} = t1
1997 val {shape
=shape2
, indexer
=indexer2
, data
=data2
} = t2
1999 if Index
.eq(shape1
,shape2
) then
2002 data
= Array
.map2 f data1 data2
}
2006 fun appi f tensor
= Array
.appi
f (toArray tensor
)
2007 fun app f tensor
= Array
.app
f (toArray tensor
)
2009 let val a
= toArray tensor
2010 in Loop
.all(0, length tensor
- 1, fn i
=>
2011 f (Array
.sub(a
, i
)))
2014 let val a
= toArray tensor
2015 in Loop
.any(0, length tensor
- 1, fn i
=>
2016 f (Array
.sub(a
, i
)))
2018 fun foldl f init tensor
= Array
.foldl f
init (toArray tensor
)
2019 fun foldln f init
{shape
, indexer
, data
=a
} index
=
2020 let val (head
,lk
,tail
) = splitList(shape
, index
)
2021 val li
= Index
.length head
2022 val lj
= Index
.length tail
2023 val c
= Array
.array(li
* lj
,init
)
2024 fun loopi (0, _
, _
) = ()
2025 |
loopi (i
, ia
, ic
) =
2026 (Array
.update(c
, ic
, f(Array
.sub(c
,ic
), Array
.sub(a
,ia
)));
2027 loopi (i
-1, ia
+1, ic
+1))
2028 fun loopk (0, ia
, _
) = ia
2029 |
loopk (k
, ia
, ic
) = (loopi (li
, ia
, ic
);
2030 loopk (k
-1, ia
+li
, ic
))
2031 fun loopj (0, _
, _
) = ()
2032 |
loopj (j
, ia
, ic
) = loopj (j
-1, loopk(lk
,ia
,ic
), ic
+li
)
2035 make
'(head @ tail
, c
)
2037 (* --- POLYMORPHIC ELEMENTWISE OPERATIONS
--- *)
2038 fun array_map
' f a
=
2039 let fun apply index
= f(Array
.sub(a
,index
)) in
2040 Tensor
.Array
.tabulate(Array
.length a
, apply
)
2042 fun map
' f t
= Tensor
.fromArray(shape t
, array_map
' f (toArray t
))
2044 let val d1
= toArray t1
2046 fun apply i
= f (Array
.sub(d1
,i
), Array
.sub(d2
,i
))
2047 val len
= Array
.length d1
2049 if Index
.eq(shape t1
, shape t2
) then
2050 Tensor
.fromArray(shape t1
, Tensor
.Array
.tabulate(len
,apply
))
2054 fun foldl
' f init
{shape
, indexer
, data
=a
} index
=
2055 let val (head
,lk
,tail
) = splitList(shape
, index
)
2056 val li
= Index
.length head
2057 val lj
= Index
.length tail
2058 val c
= Tensor
.Array
.array(li
* lj
,init
)
2059 fun loopi (0, _
, _
) = ()
2060 |
loopi (i
, ia
, ic
) =
2061 (Tensor
.Array
.update(c
,ic
,f(Tensor
.Array
.sub(c
,ic
),Array
.sub(a
,ia
)));
2062 loopi (i
-1, ia
+1, ic
+1))
2063 fun loopk (0, ia
, _
) = ia
2064 |
loopk (k
, ia
, ic
) = (loopi (li
, ia
, ic
);
2065 loopk (k
-1, ia
+li
, ic
))
2066 fun loopj (0, _
, _
) = ()
2067 |
loopj (j
, ia
, ic
) = loopj (j
-1, loopk(lk
,ia
,ic
), ic
+li
)
2070 make
'(head @ tail
, c
)
2073 end (* MonoTensor
*)
2077 LEFT INDEX CONTRACTION
:
2080 c
= c(i2
,...,in,j2
,...,jn
)
2081 = sum(a(k
,i2
,...,jn
)*b(k
,j2
,...jn
)) forall k
2082 MEANINGFUL VARIABLES
:
2087 fun do_fold_first a b c lk lj li
=
2088 let fun loopk (0, _
, _
, accum
) = accum
2089 |
loopk (k
, ia
, ib
, accum
) =
2090 let val delta
= Number
.*(Array
.sub(a
,ia
),Array
.sub(b
,ib
))
2091 in loopk (k
-1, ia
+1, ib
+1, Number
.+(delta
,accum
))
2093 fun loopj (0, ib
, ic
) = c
2094 |
loopj (j
, ib
, ic
) =
2095 let fun loopi (0, ia
, ic
) = ic
2096 |
loopi (i
, ia
, ic
) =
2097 (Array
.update(c
, ic
, loopk(lk
, ia
, ib
, Number
.zero
));
2098 loopi(i
-1, ia
+lk
, ic
+1))
2100 loopj(j
-1, ib
+lk
, loopi(li
, 0, ic
))
2106 let val (rank_a
,lk
::rest_a
,a
) = (rank ta
, shape ta
, toArray ta
)
2107 val (rank_b
,lk2
::rest_b
,b
) = (rank tb
, shape tb
, toArray tb
)
2110 else let val li
= Index
.length rest_a
2111 val lj
= Index
.length rest_b
2112 val c
= Array
.array(li
*lj
,Number
.zero
)
2113 in fromArray(rest_a @ rest_b
,
2114 do_fold_first a b c lk li lj
)
2120 LAST INDEX CONTRACTION
:
2123 c
= c(i2
,...,in,j2
,...,jn
)
2124 = sum(mult(a(i1
,i2
,...,k
),b(j1
,j2
,...,k
))) forall k
2125 MEANINGFUL VARIABLES
:
2130 fun do_fold_last a b c lk lj li
=
2131 let fun loopi (0, ia
, ic
, fac
) = ()
2132 |
loopi (i
, ia
, ic
, fac
) =
2133 let val old
= Array
.sub(c
,ic
)
2134 val inc
= Number
.*(Array
.sub(a
,ia
),fac
)
2136 Array
.update(c
,ic
,Number
.+(old
,inc
));
2137 loopi(i
-1, ia
+1, ic
+1, fac
)
2139 fun loopj (j
, ib
, ic
) =
2140 let fun loopk (0, ia
, ib
) = ()
2141 |
loopk (k
, ia
, ib
) =
2142 (loopi(li
, ia
, ic
, Array
.sub(b
,ib
));
2143 loopk(k
-1, ia
+li
, ib
+lj
))
2146 | _
=> (loopk(lk
, 0, ib
);
2147 loopj(j
-1, ib
+1, ic
+li
))
2154 let val (rank_a
,shape_a
,a
) = (rank ta
, shape ta
, toArray ta
)
2155 val (rank_b
,shape_b
,b
) = (rank tb
, shape tb
, toArray tb
)
2156 val (lk
::rest_a
) = List.rev shape_a
2157 val (lk2
::rest_b
) = List.rev shape_b
2160 else let val li
= Index
.length rest_a
2161 val lj
= Index
.length rest_b
2162 val c
= Array
.array(li
*lj
,Number
.zero
)
2163 in fromArray(List.rev rest_a @
List.rev rest_b
,
2164 do_fold_last a b c lk li lj
)
2168 (* ALGEBRAIC OPERATIONS
*)
2172 fun a
+ b
= map2 Number
.+ a b
2173 fun a
- b
= map2 Number
.- a b
2174 fun a
* b
= map2 Number
.* a b
2175 fun a
** i
= map (fn x
=> (Number
.**(x
,i
))) a
2176 fun ~ a
= map Number
.~ a
2177 fun abs a
= map Number
.abs a
2178 fun signum a
= map Number
.signum a
2179 fun a
== b
= map2
' Number
.== a b
2180 fun a
!= b
= map2
' Number
.!= a b
2181 fun toString a
= raise Domain
2182 fun fromInt a
= new([1], Number
.fromInt a
)
2183 (* TENSOR SPECIFIC OPERATIONS
*)
2184 fun *> n
= map (fn x
=> Number
.*(n
,x
))
2186 (PrettyPrint
.intList (shape t
);
2188 PrettyPrint
.sequence (length t
) appi Number
.toString t
)
2190 let fun accum (y
,x
) = Number
.max(x
,Number
.abs y
)
2191 in foldl accum Number
.zero a
2193 end (* NumberTensor
*)
2196 structure Number
= RNumber
2197 structure Array
= RNumberArray
2199 Copyright (c
) Juan Jose Garcia Ripoll
.
2200 All rights reserved
.
2201 Refer to the COPYRIGHT file for license conditions
2203 structure MonoTensor
=
2206 structure Array
= Array
2208 structure Index
= Index
2209 type elem
= Array
.elem
2210 type index
= Index
.t
2211 type tensor
= {shape
: index
, indexer
: Index
.indexer
, data
: Array
.array
}
2217 (*----- LOCALS
-----*)
2218 fun make
' (shape
, data
) =
2219 {shape
= shape
, indexer
= Index
.indexer shape
, data
= data
}
2220 fun toInt
{shape
, indexer
, data
} index
= indexer index
2221 fun splitList (l
as (a
::rest
), place
) =
2222 let fun loop (left
,here
,right
) 0 = (List.rev left
,here
,right
)
2223 |
loop (_
,_
,[]) place
= raise Index
2224 |
loop (left
,here
,a
::right
) place
=
2225 loop (here
::left
,a
,right
) (place
-1)
2228 loop ([],a
,rest
) (List.length rest
- place
)
2230 loop ([],a
,rest
) (place
- 1)
2233 (*----- STRUCTURAL OPERATIONS
& QUERIES
------*)
2234 fun new (shape
, init
) =
2235 if not (Index
.validShape shape
) then
2238 let val length
= Index
.length shape
in
2240 indexer
= Index
.indexer shape
,
2241 data
= Array
.array(length
,init
)}
2243 fun toArray
{shape
, indexer
, data
} = data
2244 fun length
{shape
, indexer
, data
} = Array
.length data
2245 fun shape
{shape
, indexer
, data
} = shape
2246 fun rank t
= List.length (shape t
)
2247 fun reshape new_shape tensor
=
2248 if Index
.validShape new_shape
then
2249 case (Index
.length new_shape
) = length tensor
of
2250 true => make
'(new_shape
, toArray tensor
)
2251 |
false => raise Match
2254 fun fromArray (s
, a
) =
2255 case Index
.validShape s
andalso
2256 ((Index
.length s
) = (Array
.length a
)) of
2258 |
false => raise Shape
2259 fun fromList (s
, a
) = fromArray (s
, Array
.fromList a
)
2260 fun tabulate (shape
,f
) =
2261 if Index
.validShape shape
then
2262 let val last
= Index
.last shape
2263 val length
= Index
.length shape
2264 val c
= Array
.array(length
, f last
)
2265 fun dotable (c
, indices
, i
) =
2266 (Array
.update(c
, i
, f indices
);
2269 else dotable(c
, Index
.prev
' shape indices
, i
-1))
2270 in make
'(shape
,dotable(c
, Index
.prev
' shape last
, length
-2))
2274 (*----- ELEMENTWISE OPERATIONS
-----*)
2275 fun sub (t
, index
) = Array
.sub(#data t
, toInt t index
)
2276 fun update (t
, index
, value
) =
2277 Array
.update(toArray t
, toInt t index
, value
)
2278 fun map f
{shape
, indexer
, data
} =
2279 {shape
= shape
, indexer
= indexer
, data
= Array
.map f data
}
2281 let val {shape
=shape1
, indexer
=indexer1
, data
=data1
} = t1
2282 val {shape
=shape2
, indexer
=indexer2
, data
=data2
} = t2
2284 if Index
.eq(shape1
,shape2
) then
2287 data
= Array
.map2 f data1 data2
}
2291 fun appi f tensor
= Array
.appi
f (toArray tensor
)
2292 fun app f tensor
= Array
.app
f (toArray tensor
)
2294 let val a
= toArray tensor
2295 in Loop
.all(0, length tensor
- 1, fn i
=>
2296 f (Array
.sub(a
, i
)))
2299 let val a
= toArray tensor
2300 in Loop
.any(0, length tensor
- 1, fn i
=>
2301 f (Array
.sub(a
, i
)))
2303 fun foldl f init tensor
= Array
.foldl f
init (toArray tensor
)
2304 fun foldln f init
{shape
, indexer
, data
=a
} index
=
2305 let val (head
,lk
,tail
) = splitList(shape
, index
)
2306 val li
= Index
.length head
2307 val lj
= Index
.length tail
2308 val c
= Array
.array(li
* lj
,init
)
2309 fun loopi (0, _
, _
) = ()
2310 |
loopi (i
, ia
, ic
) =
2311 (Array
.update(c
, ic
, f(Array
.sub(c
,ic
), Array
.sub(a
,ia
)));
2312 loopi (i
-1, ia
+1, ic
+1))
2313 fun loopk (0, ia
, _
) = ia
2314 |
loopk (k
, ia
, ic
) = (loopi (li
, ia
, ic
);
2315 loopk (k
-1, ia
+li
, ic
))
2316 fun loopj (0, _
, _
) = ()
2317 |
loopj (j
, ia
, ic
) = loopj (j
-1, loopk(lk
,ia
,ic
), ic
+li
)
2320 make
'(head @ tail
, c
)
2322 (* --- POLYMORPHIC ELEMENTWISE OPERATIONS
--- *)
2323 fun array_map
' f a
=
2324 let fun apply index
= f(Array
.sub(a
,index
)) in
2325 Tensor
.Array
.tabulate(Array
.length a
, apply
)
2327 fun map
' f t
= Tensor
.fromArray(shape t
, array_map
' f (toArray t
))
2329 let val d1
= toArray t1
2331 fun apply i
= f (Array
.sub(d1
,i
), Array
.sub(d2
,i
))
2332 val len
= Array
.length d1
2334 if Index
.eq(shape t1
, shape t2
) then
2335 Tensor
.fromArray(shape t1
, Tensor
.Array
.tabulate(len
,apply
))
2339 fun foldl
' f init
{shape
, indexer
, data
=a
} index
=
2340 let val (head
,lk
,tail
) = splitList(shape
, index
)
2341 val li
= Index
.length head
2342 val lj
= Index
.length tail
2343 val c
= Tensor
.Array
.array(li
* lj
,init
)
2344 fun loopi (0, _
, _
) = ()
2345 |
loopi (i
, ia
, ic
) =
2346 (Tensor
.Array
.update(c
,ic
,f(Tensor
.Array
.sub(c
,ic
),Array
.sub(a
,ia
)));
2347 loopi (i
-1, ia
+1, ic
+1))
2348 fun loopk (0, ia
, _
) = ia
2349 |
loopk (k
, ia
, ic
) = (loopi (li
, ia
, ic
);
2350 loopk (k
-1, ia
+li
, ic
))
2351 fun loopj (0, _
, _
) = ()
2352 |
loopj (j
, ia
, ic
) = loopj (j
-1, loopk(lk
,ia
,ic
), ic
+li
)
2355 make
'(head @ tail
, c
)
2358 end (* MonoTensor
*)
2362 LEFT INDEX CONTRACTION
:
2365 c
= c(i2
,...,in,j2
,...,jn
)
2366 = sum(a(k
,i2
,...,jn
)*b(k
,j2
,...jn
)) forall k
2367 MEANINGFUL VARIABLES
:
2372 fun do_fold_first a b c lk lj li
=
2373 let fun loopk (0, _
, _
, accum
) = accum
2374 |
loopk (k
, ia
, ib
, accum
) =
2375 let val delta
= Number
.*(Array
.sub(a
,ia
),Array
.sub(b
,ib
))
2376 in loopk (k
-1, ia
+1, ib
+1, Number
.+(delta
,accum
))
2378 fun loopj (0, ib
, ic
) = c
2379 |
loopj (j
, ib
, ic
) =
2380 let fun loopi (0, ia
, ic
) = ic
2381 |
loopi (i
, ia
, ic
) =
2382 (Array
.update(c
, ic
, loopk(lk
, ia
, ib
, Number
.zero
));
2383 loopi(i
-1, ia
+lk
, ic
+1))
2385 loopj(j
-1, ib
+lk
, loopi(li
, 0, ic
))
2391 let val (rank_a
,lk
::rest_a
,a
) = (rank ta
, shape ta
, toArray ta
)
2392 val (rank_b
,lk2
::rest_b
,b
) = (rank tb
, shape tb
, toArray tb
)
2395 else let val li
= Index
.length rest_a
2396 val lj
= Index
.length rest_b
2397 val c
= Array
.array(li
*lj
,Number
.zero
)
2398 in fromArray(rest_a @ rest_b
,
2399 do_fold_first a b c lk li lj
)
2405 LAST INDEX CONTRACTION
:
2408 c
= c(i2
,...,in,j2
,...,jn
)
2409 = sum(mult(a(i1
,i2
,...,k
),b(j1
,j2
,...,k
))) forall k
2410 MEANINGFUL VARIABLES
:
2415 fun do_fold_last a b c lk lj li
=
2416 let fun loopi (0, ia
, ic
, fac
) = ()
2417 |
loopi (i
, ia
, ic
, fac
) =
2418 let val old
= Array
.sub(c
,ic
)
2419 val inc
= Number
.*(Array
.sub(a
,ia
),fac
)
2421 Array
.update(c
,ic
,Number
.+(old
,inc
));
2422 loopi(i
-1, ia
+1, ic
+1, fac
)
2424 fun loopj (j
, ib
, ic
) =
2425 let fun loopk (0, ia
, ib
) = ()
2426 |
loopk (k
, ia
, ib
) =
2427 (loopi(li
, ia
, ic
, Array
.sub(b
,ib
));
2428 loopk(k
-1, ia
+li
, ib
+lj
))
2431 | _
=> (loopk(lk
, 0, ib
);
2432 loopj(j
-1, ib
+1, ic
+li
))
2439 let val (rank_a
,shape_a
,a
) = (rank ta
, shape ta
, toArray ta
)
2440 val (rank_b
,shape_b
,b
) = (rank tb
, shape tb
, toArray tb
)
2441 val (lk
::rest_a
) = List.rev shape_a
2442 val (lk2
::rest_b
) = List.rev shape_b
2445 else let val li
= Index
.length rest_a
2446 val lj
= Index
.length rest_b
2447 val c
= Array
.array(li
*lj
,Number
.zero
)
2448 in fromArray(List.rev rest_a @
List.rev rest_b
,
2449 do_fold_last a b c lk li lj
)
2453 (* ALGEBRAIC OPERATIONS
*)
2457 fun a
+ b
= map2 Number
.+ a b
2458 fun a
- b
= map2 Number
.- a b
2459 fun a
* b
= map2 Number
.* a b
2460 fun a
** i
= map (fn x
=> (Number
.**(x
,i
))) a
2461 fun ~ a
= map Number
.~ a
2462 fun abs a
= map Number
.abs a
2463 fun signum a
= map Number
.signum a
2464 fun a
== b
= map2
' Number
.== a b
2465 fun a
!= b
= map2
' Number
.!= a b
2466 fun toString a
= raise Domain
2467 fun fromInt a
= new([1], Number
.fromInt a
)
2468 (* TENSOR SPECIFIC OPERATIONS
*)
2469 fun *> n
= map (fn x
=> Number
.*(n
,x
))
2471 (PrettyPrint
.intList (shape t
);
2473 PrettyPrint
.sequence (length t
) appi Number
.toString t
)
2474 fun a
/ b
= map2 Number
./ a b
2475 fun recip a
= map Number
.recip a
2476 fun ln a
= map Number
.ln a
2477 fun pow (a
, b
) = map (fn x
=> (Number
.pow(x
,b
))) a
2478 fun exp a
= map Number
.exp a
2479 fun sqrt a
= map Number
.sqrt a
2480 fun cos a
= map Number
.cos a
2481 fun sin a
= map Number
.sin a
2482 fun tan a
= map Number
.tan a
2483 fun sinh a
= map Number
.sinh a
2484 fun cosh a
= map Number
.cosh a
2485 fun tanh a
= map Number
.tanh a
2486 fun asin a
= map Number
.asin a
2487 fun acos a
= map Number
.acos a
2488 fun atan a
= map Number
.atan a
2489 fun asinh a
= map Number
.asinh a
2490 fun acosh a
= map Number
.acosh a
2491 fun atanh a
= map Number
.atanh a
2492 fun atan2 (a
,b
) = map2 Number
.atan2 a b
2494 let fun accum (y
,x
) = Number
.max(x
,Number
.abs y
)
2495 in foldl accum Number
.zero a
2498 let fun accum (y
,x
) = Number
.+(x
,Number
.abs y
)
2499 in foldl accum Number
.zero a
2502 let fun accum (y
,x
) = Number
.+(x
, Number
.*(y
,y
))
2503 in Number
.sqrt(foldl accum Number
.zero a
)
2508 structure Number
= CNumber
2509 structure Array
= CNumberArray
2511 Copyright (c
) Juan Jose Garcia Ripoll
.
2512 All rights reserved
.
2513 Refer to the COPYRIGHT file for license conditions
2515 structure MonoTensor
=
2518 structure Array
= Array
2520 structure Index
= Index
2521 type elem
= Array
.elem
2522 type index
= Index
.t
2523 type tensor
= {shape
: index
, indexer
: Index
.indexer
, data
: Array
.array
}
2529 (*----- LOCALS
-----*)
2530 fun make
' (shape
, data
) =
2531 {shape
= shape
, indexer
= Index
.indexer shape
, data
= data
}
2532 fun toInt
{shape
, indexer
, data
} index
= indexer index
2533 fun splitList (l
as (a
::rest
), place
) =
2534 let fun loop (left
,here
,right
) 0 = (List.rev left
,here
,right
)
2535 |
loop (_
,_
,[]) place
= raise Index
2536 |
loop (left
,here
,a
::right
) place
=
2537 loop (here
::left
,a
,right
) (place
-1)
2540 loop ([],a
,rest
) (List.length rest
- place
)
2542 loop ([],a
,rest
) (place
- 1)
2545 (*----- STRUCTURAL OPERATIONS
& QUERIES
------*)
2546 fun new (shape
, init
) =
2547 if not (Index
.validShape shape
) then
2550 let val length
= Index
.length shape
in
2552 indexer
= Index
.indexer shape
,
2553 data
= Array
.array(length
,init
)}
2555 fun toArray
{shape
, indexer
, data
} = data
2556 fun length
{shape
, indexer
, data
} = Array
.length data
2557 fun shape
{shape
, indexer
, data
} = shape
2558 fun rank t
= List.length (shape t
)
2559 fun reshape new_shape tensor
=
2560 if Index
.validShape new_shape
then
2561 case (Index
.length new_shape
) = length tensor
of
2562 true => make
'(new_shape
, toArray tensor
)
2563 |
false => raise Match
2566 fun fromArray (s
, a
) =
2567 case Index
.validShape s
andalso
2568 ((Index
.length s
) = (Array
.length a
)) of
2570 |
false => raise Shape
2571 fun fromList (s
, a
) = fromArray (s
, Array
.fromList a
)
2572 fun tabulate (shape
,f
) =
2573 if Index
.validShape shape
then
2574 let val last
= Index
.last shape
2575 val length
= Index
.length shape
2576 val c
= Array
.array(length
, f last
)
2577 fun dotable (c
, indices
, i
) =
2578 (Array
.update(c
, i
, f indices
);
2581 else dotable(c
, Index
.prev
' shape indices
, i
-1))
2582 in make
'(shape
,dotable(c
, Index
.prev
' shape last
, length
-2))
2586 (*----- ELEMENTWISE OPERATIONS
-----*)
2587 fun sub (t
, index
) = Array
.sub(#data t
, toInt t index
)
2588 fun update (t
, index
, value
) =
2589 Array
.update(toArray t
, toInt t index
, value
)
2590 fun map f
{shape
, indexer
, data
} =
2591 {shape
= shape
, indexer
= indexer
, data
= Array
.map f data
}
2593 let val {shape
=shape1
, indexer
=indexer1
, data
=data1
} = t1
2594 val {shape
=shape2
, indexer
=indexer2
, data
=data2
} = t2
2596 if Index
.eq(shape1
,shape2
) then
2599 data
= Array
.map2 f data1 data2
}
2603 fun appi f tensor
= Array
.appi
f (toArray tensor
, 0, NONE
)
2604 fun app f tensor
= Array
.app
f (toArray tensor
)
2606 let val a
= toArray tensor
2607 in Loop
.all(0, length tensor
- 1, fn i
=>
2608 f (Array
.sub(a
, i
)))
2611 let val a
= toArray tensor
2612 in Loop
.any(0, length tensor
- 1, fn i
=>
2613 f (Array
.sub(a
, i
)))
2615 fun foldl f init tensor
= Array
.foldl f
init (toArray tensor
)
2616 fun foldln f init
{shape
, indexer
, data
=a
} index
=
2617 let val (head
,lk
,tail
) = splitList(shape
, index
)
2618 val li
= Index
.length head
2619 val lj
= Index
.length tail
2620 val c
= Array
.array(li
* lj
,init
)
2621 fun loopi (0, _
, _
) = ()
2622 |
loopi (i
, ia
, ic
) =
2623 (Array
.update(c
, ic
, f(Array
.sub(c
,ic
), Array
.sub(a
,ia
)));
2624 loopi (i
-1, ia
+1, ic
+1))
2625 fun loopk (0, ia
, _
) = ia
2626 |
loopk (k
, ia
, ic
) = (loopi (li
, ia
, ic
);
2627 loopk (k
-1, ia
+li
, ic
))
2628 fun loopj (0, _
, _
) = ()
2629 |
loopj (j
, ia
, ic
) = loopj (j
-1, loopk(lk
,ia
,ic
), ic
+li
)
2632 make
'(head @ tail
, c
)
2634 (* --- POLYMORPHIC ELEMENTWISE OPERATIONS
--- *)
2635 fun array_map
' f a
=
2636 let fun apply index
= f(Array
.sub(a
,index
)) in
2637 Tensor
.Array
.tabulate(Array
.length a
, apply
)
2639 fun map
' f t
= Tensor
.fromArray(shape t
, array_map
' f (toArray t
))
2641 let val d1
= toArray t1
2643 fun apply i
= f (Array
.sub(d1
,i
), Array
.sub(d2
,i
))
2644 val len
= Array
.length d1
2646 if Index
.eq(shape t1
, shape t2
) then
2647 Tensor
.fromArray(shape t1
, Tensor
.Array
.tabulate(len
,apply
))
2651 fun foldl
' f init
{shape
, indexer
, data
=a
} index
=
2652 let val (head
,lk
,tail
) = splitList(shape
, index
)
2653 val li
= Index
.length head
2654 val lj
= Index
.length tail
2655 val c
= Tensor
.Array
.array(li
* lj
,init
)
2656 fun loopi (0, _
, _
) = ()
2657 |
loopi (i
, ia
, ic
) =
2658 (Tensor
.Array
.update(c
,ic
,f(Tensor
.Array
.sub(c
,ic
),Array
.sub(a
,ia
)));
2659 loopi (i
-1, ia
+1, ic
+1))
2660 fun loopk (0, ia
, _
) = ia
2661 |
loopk (k
, ia
, ic
) = (loopi (li
, ia
, ic
);
2662 loopk (k
-1, ia
+li
, ic
))
2663 fun loopj (0, _
, _
) = ()
2664 |
loopj (j
, ia
, ic
) = loopj (j
-1, loopk(lk
,ia
,ic
), ic
+li
)
2667 make
'(head @ tail
, c
)
2670 end (* MonoTensor
*)
2674 LEFT INDEX CONTRACTION
:
2677 c
= c(i2
,...,in,j2
,...,jn
)
2678 = sum(a(k
,i2
,...,jn
)*b(k
,j2
,...jn
)) forall k
2679 MEANINGFUL VARIABLES
:
2684 fun do_fold_first a b c lk lj li
=
2685 let fun loopk (0, _
, _
, r
, i
) = Number
.make(r
,i
)
2686 |
loopk (k
, ia
, ib
, r
, i
) =
2687 let val (ar
, ai
) = Array
.sub(a
,ia
)
2688 val (br
, bi
) = Array
.sub(b
,ib
)
2689 val dr
= ar
* br
- ai
* bi
2690 val di
= ar
* bi
+ ai
* br
2691 in loopk (k
-1, ia
+1, ib
+1, r
+dr
, i
+di
)
2693 fun loopj (0, ib
, ic
) = c
2694 |
loopj (j
, ib
, ic
) =
2695 let fun loopi (0, ia
, ic
) = ic
2696 |
loopi (i
, ia
, ic
) =
2697 (Array
.update(c
, ic
, loopk(lk
, ia
, ib
, RNumber
.zero
, RNumber
.zero
));
2698 loopi(i
-1, ia
+lk
, ic
+1))
2699 in loopj(j
-1, ib
+lk
, loopi(li
, 0, ic
))
2705 let val (rank_a
,lk
::rest_a
,a
) = (rank ta
, shape ta
, toArray ta
)
2706 val (rank_b
,lk2
::rest_b
,b
) = (rank tb
, shape tb
, toArray tb
)
2709 else let val li
= Index
.length rest_a
2710 val lj
= Index
.length rest_b
2711 val c
= Array
.array(li
*lj
,Number
.zero
)
2712 in fromArray(rest_a @ rest_b
, do_fold_first a b c lk li lj
)
2718 LAST INDEX CONTRACTION
:
2721 c
= c(i2
,...,in,j2
,...,jn
)
2722 = sum(mult(a(i1
,i2
,...,k
),b(j1
,j2
,...,k
))) forall k
2723 MEANINGFUL VARIABLES
:
2728 fun do_fold_last a b c lk lj li
=
2729 let fun loopi(0, _
, _
, _
, _
) = ()
2730 |
loopi(i
, ia
, ic
, br
, bi
) =
2731 let val (cr
,ci
) = Array
.sub(c
,ic
)
2732 val (ar
,ai
) = Array
.sub(a
,ia
)
2733 val dr
= (ar
* br
- ai
* bi
)
2734 val di
= (ar
* bi
+ ai
* br
)
2736 Array
.update(c
,ic
,Number
.make(cr
+dr
,ci
+di
));
2737 loopi(i
-1, ia
+1, ic
+1, br
, bi
)
2739 fun loopj(j
, ib
, ic
) =
2740 let fun loopk(0, _
, _
) = ()
2741 |
loopk(k
, ia
, ib
) =
2742 let val (br
, bi
) = Array
.sub(b
,ib
)
2744 loopi(li
, ia
, ic
, br
, bi
);
2745 loopk(k
-1, ia
+li
, ib
+lj
)
2749 | _
=> (loopk(lk
, 0, ib
);
2750 loopj(j
-1, ib
+1, ic
+li
))
2757 let val (rank_a
,shape_a
,a
) = (rank ta
, shape ta
, toArray ta
)
2758 val (rank_b
,shape_b
,b
) = (rank tb
, shape tb
, toArray tb
)
2759 val (lk
::rest_a
) = List.rev shape_a
2760 val (lk2
::rest_b
) = List.rev shape_b
2762 if not(lk
= lk2
) then
2765 let val li
= Index
.length rest_a
2766 val lj
= Index
.length rest_b
2767 val c
= Array
.array(li
*lj
,Number
.zero
)
2769 fromArray(List.rev rest_a @
List.rev rest_b
,
2770 do_fold_last a b c lk li lj
)
2774 (* ALGEBRAIC OPERATIONS
*)
2778 fun a
+ b
= map2 Number
.+ a b
2779 fun a
- b
= map2 Number
.- a b
2780 fun a
* b
= map2 Number
.* a b
2781 fun a
** i
= map (fn x
=> (Number
.**(x
,i
))) a
2782 fun ~ a
= map Number
.~ a
2783 fun abs a
= map Number
.abs a
2784 fun signum a
= map Number
.signum a
2785 fun a
== b
= map2
' Number
.== a b
2786 fun a
!= b
= map2
' Number
.!= a b
2787 fun toString a
= raise Domain
2788 fun fromInt a
= new([1], Number
.fromInt a
)
2789 (* TENSOR SPECIFIC OPERATIONS
*)
2790 fun *> n
= map (fn x
=> Number
.*(n
,x
))
2792 (PrettyPrint
.intList (shape t
);
2794 PrettyPrint
.sequence (length t
) appi Number
.toString t
)
2795 fun a
/ b
= map2 Number
./ a b
2796 fun recip a
= map Number
.recip a
2797 fun ln a
= map Number
.ln a
2798 fun pow (a
, b
) = map (fn x
=> (Number
.pow(x
,b
))) a
2799 fun exp a
= map Number
.exp a
2800 fun sqrt a
= map Number
.sqrt a
2801 fun cos a
= map Number
.cos a
2802 fun sin a
= map Number
.sin a
2803 fun tan a
= map Number
.tan a
2804 fun sinh a
= map Number
.sinh a
2805 fun cosh a
= map Number
.cosh a
2806 fun tanh a
= map Number
.tanh a
2807 fun asin a
= map Number
.asin a
2808 fun acos a
= map Number
.acos a
2809 fun atan a
= map Number
.atan a
2810 fun asinh a
= map Number
.asinh a
2811 fun acosh a
= map Number
.acosh a
2812 fun atanh a
= map Number
.atanh a
2813 fun atan2 (a
,b
) = map2 Number
.atan2 a b
2815 let fun accum (y
,x
) = RNumber
.max(x
, Number
.realPart(Number
.abs y
))
2816 in foldl accum RNumber
.zero a
2819 let fun accum (y
,x
) = RNumber
.+(x
, Number
.realPart(Number
.abs y
))
2820 in foldl accum RNumber
.zero a
2823 let fun accum (y
,x
) = RNumber
.+(x
, Number
.abs2 y
)
2824 in RNumber
.sqrt(foldl accum RNumber
.zero a
)
2827 structure MathFile
=
2830 type file
= TextIO.instream
2834 fun assert NONE
= raise Data
2835 |
assert (SOME a
) = a
2837 (* ------------------ INPUT
--------------------- *)
2839 fun intRead file
= assert(TextIO.scanStream INumber
.scan file
)
2840 fun realRead file
= assert(TextIO.scanStream RNumber
.scan file
)
2841 fun complexRead file
= assert(TextIO.scanStream CNumber
.scan file
)
2843 fun listRead eltScan file
=
2844 let val length
= intRead file
2845 fun eltRead file
= assert(TextIO.scanStream eltScan file
)
2846 fun loop (0,accum
) = accum
2847 |
loop (i
,accum
) = loop(i
-1, eltRead file
:: accum
)
2851 else List.rev(loop(length
,[]))
2854 fun intListRead file
= listRead INumber
.scan file
2855 fun realListRead file
= listRead RNumber
.scan file
2856 fun complexListRead file
= listRead CNumber
.scan file
2858 fun intTensorRead file
=
2859 let val shape
= intListRead file
2860 val length
= Index
.length shape
2861 val first
= intRead file
2862 val a
= ITensor
.Array
.array(length
, first
)
2863 fun loop
0 = ITensor
.fromArray(shape
, a
)
2864 | loop j
= (ITensor
.Array
.update(a
, length
-j
, intRead file
);
2866 in loop (length
- 1)
2869 fun realTensorRead file
=
2870 let val shape
= intListRead file
2871 val length
= Index
.length shape
2872 val first
= realRead file
2873 val a
= RTensor
.Array
.array(length
, first
)
2874 fun loop
0 = RTensor
.fromArray(shape
, a
)
2875 | loop j
= (RTensor
.Array
.update(a
, length
-j
, realRead file
);
2877 in loop (length
- 1)
2880 fun complexTensorRead file
=
2881 let val shape
= intListRead file
2882 val length
= Index
.length shape
2883 val first
= complexRead file
2884 val a
= CTensor
.Array
.array(length
, first
)
2885 fun loop j
= if j
= length
2886 then CTensor
.fromArray(shape
, a
)
2887 else (CTensor
.Array
.update(a
, j
, complexRead file
);
2892 (* ------------------ OUTPUT
-------------------- *)
2893 fun linedOutput(file
, x
) = (TextIO.output(file
, x
); TextIO.output(file
, "\n"))
2895 fun intWrite file x
= linedOutput(file
, INumber
.toString x
)
2896 fun realWrite file x
= linedOutput(file
, RNumber
.toString x
)
2897 fun complexWrite file x
=
2898 let val (r
,i
) = CNumber
.split x
2899 in linedOutput(file
, concat
[RNumber
.toString r
, " ", RNumber
.toString i
])
2902 fun listWrite converter file x
=
2903 (intWrite
file (length x
);
2904 List.app (fn x
=> (linedOutput(file
, converter x
))) x
)
2906 fun intListWrite file x
= listWrite INumber
.toString file x
2907 fun realListWrite file x
= listWrite RNumber
.toString file x
2908 fun complexListWrite file x
= listWrite CNumber
.toString file x
2910 fun intTensorWrite file x
= (intListWrite
file (ITensor
.shape x
); ITensor
.app (fn x
=> (intWrite file x
)) x
)
2911 fun realTensorWrite file x
= (intListWrite
file (RTensor
.shape x
); RTensor
.app (fn x
=> (realWrite file x
)) x
)
2912 fun complexTensorWrite file x
= (intListWrite
file (CTensor
.shape x
); CTensor
.app (fn x
=> (complexWrite file x
)) x
)
2916 | loop n f
= (f(); loop (n
-1) f
)
2918 fun test_operator new list_op list_sizes
=
2919 let fun test_many list_op size
=
2920 let fun test_op (times
,f
) =
2921 let val a
= new size
2922 in (EvalTimer
.timerOn();
2923 loop
times (fn _
=> f(a
,a
));
2924 let val t
= LargeInt
.toInt(EvalTimer
.timerRead()) div times
2925 val i
= StringCvt.padLeft #
" " 6 (Int.toString t
)
2930 print (Int.toString size
);
2932 List.app test_op list_op
;
2935 in List.app (test_many list_op
) list_sizes
2943 let val operators
= [(20, RTensor
.+), (20, RTensor
.* ), (20, RTensor
./),
2944 (4, fn (a
,b
) => RTensor
.+* a b
),
2945 (4, fn (a
,b
) => RTensor
.*+ a b
)]
2946 fun constructor size
= RTensor
.new([size
,size
],1.0)
2948 print
"Real tensors: (+, *, /, +*, *+)\n";
2949 test_operator constructor operators
[100,200,300,400,500];
2954 let val operators
= [(20, CTensor
.+), (20, CTensor
.* ), (20, CTensor
./),
2955 (4, fn (a
,b
) => CTensor
.+* a b
),
2956 (4, fn (a
,b
) => CTensor
.*+ a b
)]
2957 fun constructor size
= CTensor
.new([size
,size
],CNumber
.one
)
2959 print
"Real tensors: (+, *, /, +*, *+)\n";
2960 test_operator constructor operators
[100,200,300,400,500];