Home My Page Projects Code Snippets Project Openings SML/NJ
Summary Activity Forums Tracker Lists Tasks Docs Surveys News SCM Files

SCM Repository

[smlnj] Annotation of /smlnj-lib/trunk/Util/real-order-stats.sml
ViewVC logotype

Annotation of /smlnj-lib/trunk/Util/real-order-stats.sml

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1720 - (view) (download)
Original Path: sml/trunk/src/smlnj-lib/Util/real-order-stats.sml

1 : mblume 1720 (* real-order-stats.sml
2 :     *
3 :     * Randomized linear-time selection from an unordered sample.
4 :     *
5 :     * Copyright (c) 2004 by The Fellowship of SML/NJ
6 :     *
7 :     * Author: Matthias Blume (blume@tti-c.org)
8 :     *)
9 :     structure RealOrderStats : sig
10 :    
11 :     (* select the i-th order statistic *)
12 :     val select : real array * int -> real
13 :     val select' : real ArraySlice.slice * int -> real
14 :    
15 :     (* calculate the median:
16 :     * if N is odd, then this is the (floor(N/2))th order statistic
17 :     * otherwise it is the average of (N/2-1)th and (N/2)th *)
18 :     val median : real array -> real
19 :     val median' : real ArraySlice.slice -> real
20 :    
21 :     end = struct
22 :    
23 :     infix 8 $ val op $ = Unsafe.Array.sub
24 :     infix 3 <- fun (a, i) <- x = Unsafe.Array.update (a, i, x)
25 :    
26 :     (* initialize random number generator *)
27 :     val rand = Random.rand (123, 73256)
28 :    
29 :     (* select i-th order statistic from unsorted array with
30 :     * starting point p and ending point r (inclusive): *)
31 :     fun select0 (a: real array, p, r, i) =
32 :     let fun x + y = Word.toIntX (Word.+ (Word.fromInt x, Word.fromInt y))
33 :     fun x - y = Word.toIntX (Word.- (Word.fromInt x, Word.fromInt y))
34 :     (* random partition: *)
35 :     fun rp (p, r) =
36 :     let fun sw(i,j) = let val t=a$i in (a,i)<-a$j; (a,j)<-t end
37 :     val q = Random.randRange (p, r) rand
38 :     val qv = a$q
39 :     val _ = if q<>p then ((a,q)<-a$p; (a,p)<-qv) else ()
40 :     fun up i = if i>r orelse qv < a$i then i else up(i+1)
41 :     fun dn i = if i>=p andalso qv < a$i then dn(i-1) else i
42 :     fun lp (i, j) =
43 :     let val (i, j) = (up i, dn j)
44 :     in if i>j then let val q' = i-1 in sw(p,q'); (q',qv) end
45 :     else (sw(i,j); lp (i+1, j-1))
46 :     end
47 :     in lp (p+1, r) end
48 :     (* random select: *)
49 :     fun rs (p, r) =
50 :     if p=r then a$r
51 :     else let val (q, qv) = rp (p, r)
52 :     in if i=q then qv else if i<q then rs(p,q-1) else rs(q+1,r)
53 :     end
54 :     in rs (p, r) end
55 :    
56 :     fun select (a, i) = select0 (a, 0, Array.length a - 1, i)
57 :     fun select' (s, i) = let val (a, p, l) = ArraySlice.base s
58 :     in select0 (a, p, p+l-1, p+i) end
59 :    
60 :     fun median0 (a, p, len) =
61 :     let val mid = p + len div 2
62 :     val r = p + len - 1
63 :     val m0 = select0 (a, p, r, mid)
64 :     in if len mod 2 = 1 then m0
65 :     else let fun l(i,m) = if i>=mid then m else l(i+1, Real.max(a$i,m))
66 :     in (l(p+1,a$p) + m0) / 2.0 end
67 :     end
68 :    
69 :     fun median a = median0 (a, 0, Array.length a)
70 :     fun median' s = median0 (ArraySlice.base s)
71 :     end

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0