summaryrefslogtreecommitdiff
path: root/web/funnelAC/contrib/dek.fw
blob: a034f137ec080f2bdb64f823239f68bd03184cd2 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
@p typesetter = tex

% Set your own page size here . . . .

\hsize 160mm
\vsize 245mm

\def\fwseca#1#2{\fwlibc{#1}{#2}}
\def\fwsecb#1#2{\fwlibd{#1}{#2}}
\def\fwsecc#1#2{\fwlibe{#1}{#2}}
\def\fwsecd#1#2{\fwlibe{#1}{#2}}
\def\fwsece#1#2{\fwlibe{#1}{#2}}

@t vskip 40 mm
@t title titlefont centre "A Programming Pearl"
@t title titlefont centre "Generating sorted random numbers"
@t vskip 20 mm
@t title smalltitlefont centre "Don Knuth"
@t vskip 10 mm
@t title normalfont centre "FunnelWeb by Richard Walker"
@t new_page
@t table_of_contents
@t new_page

@!------------------------------------------------------------
@a@<Introduction@>

Jon Bentley recently discussed the following interesting problem
as one of his `Programming Pearls' [{\sl Communications of the
ACM 27\/} (December, 1984), 1179--1182]:

The input consists of two integers @{M@} and @{N@}, with @{M@}$<$@{N@}.
The output is a sorted list of @{M@} random numbers in the range
1..@{N@} in which no integer occurs more than once.  For probability
buffs, we desire a sorted selection without replacement in which each
selection occurs equiprobably.

The present program illustrates what I think is the best solution to
the problem, when @{M@} is reasonably large yet small compared to @{N@}.
It's the method described tersely in the answer to exercise 3.4.2--15 of
my book {\sl Seminumerical Algorithms}, pp.~141 and 555.

@!------------------------------------------------------------
@a@<Some macros@>

For simplicity, all input and output in this program is assumed to be
handled at the terminal.  The @{WEB@} macros @{read_terminal@},
@{print@}, and @{print_ln@} defined here can easily be changed to
accommodate other conventions.

Input a value from the terminal.
@$@<read_terminal@>@(@1@)@M==@{read(@1)@}

Output to the terminal.
@$@<print@>@(@1@)@M==@{write(@1)@}

Output to the terminal and end the line.
@$@<print_ln@>@(@1@)@M==@{writeln(@1)@}

@!------------------------------------------------------------
@a@<An outline of the program@>

Here's an outline of the entire Pascal program:

@o@<sample.p@>==@{@-
program sample(input,output);
  var @<Global variables@>
  @<The random number generation procedure@>
  begin @<The main program@>;
  end.
@}

@!------------------------------------------------------------
@a@<Global variables@>

The global variables @{M@} and @{N@} have already been mentioned; we had
better declare them.  Other global variables will be declared later.

@{M_max@} is the maximum value of @{M@} allowed in this program.
@$@<M_max@>@M==@{5000@}

@{2M-1_max@} is $2M-1$ for use later on.
@$@<2M-1_max@>@M==@{9999@}

@$@<Global variables@>+=@{@-
M: integer; { size of the sample }
N: integer; { size of the population }
@}

@!------------------------------------------------------------
@a@<Random number generator@>

We assume the existence of a system routine called @{rand_int(i,j)@}
that returns a random integer chosen uniformly in the range $i..j$.

@$@<The random number generation procedure@>==@{@-
function rand_int(i,j: integer): integer; external;
@}

@!------------------------------------------------------------
@a@<A plan of attack@>

After the user has specified @{M@} and @{N@}, we compute the sample by
following a general procedure recommended by Bentley:

@$@<The main program@>==@{@-
@<Establish the values of M and N@>;
size := 0; @<Initialise set S to empty@>;
while size < M do
  begin
    T := rand_int(1,N);
    @<If T is not in S, insert it and increase size@>;
  end;
@<Print the elements of S in sorted order@>@}

@!------------------------------------------------------------
@a@<More globals@>

The main program just sketched has introduced several more
globals.  There's a set @{S@} of integers, whose representation
will be deferred until later; but we can declare two auxiliary
integer variables now.

@$@<Global variables@>+=@{@-
size: integer; { the number of elements in set S }
T: integer; { new candidate for membership in S }
@}

@!------------------------------------------------------------
@a

The first order of business is to have a short dialogue with the
user.

@$@<Establish the values of M and N@>==@{@-
repeat @<print@>@('population size: N = '@);
  @<read_terminal@>@(N@);
  if N <= 0 then
    @<print_ln@>@('N should be positive!'@);
until N > 0;
repeat @<print@>@('sample size: M = '@);
  @<read_terminal@>@(M@);
  if M < 0 then
    @<print_ln@>@('M shouldn''t be negative!'@)
  else if M > N then
    @<print_ln@>@('M shouldn''t exceed N!'@)
  else if M > @<M_max@> then
    @<print_ln@>@('(Sorry, M must be at most ',@<M_max@>:1,'.)'@);
until (M >= 0) and (M <= N) and (M <= @<M_max@>)@}

@!------------------------------------------------------------
@a@<An ordered hash table@>

The key idea to an efficient solution of this sampling problem is
to maintain a set whose entries are easily sorted.  The method of
`ordered hash tables' [Amble and Knuth, {\sl The Computer Journal
17\/} (May 1974), 135--142] is ideally suited to this task, as we
shall see.

Ordered hashing is similar to ordinary linear probing, except that
the relative order of keys is taken into account.  The cited
paper derives theoretical results that will not be rederived
here, but we shall use the following fundamental property:  {\sl
The entries of an ordered hash table are independent of the order
in which its keys were inserted}.  Thus, an ordered hash table is
a `canonical' representation of its set of entries.

We shall represent @{S@} by an array of $2M$ integers.  Since
Pascal doesn't permit arrays of variable size, we must leave room
for the largest possible table.

@$@<Global variables@>+=@{@-
hash: array[0..@<2M-1_max@>] of integer;
    { the ordered hash table }
H: 0..@<2M-1_max@>; { an index into hash }
H_max: 0..@<2M-1_max@>; { the current hash size }
alpha: real; { the ratio of table size to N }
@}

@!------------------------------------------------------------
@a

@$@<Initialise set S to empty@>==@{@-
H_max := 2 * M - 1; alpha := 2 * M / N;
for H := 0 to H_max do hash[H] := 0@}

@!------------------------------------------------------------
@a

Now we come to the interesting part, where the algorithm tries to
insert @{T@} into an ordered hash table.  The hash address
$H=\lfloor2M(T-1)/N\rfloor$ is used as a starting point, since
this quantity is monotonic in @{T@} and almost uniformly
distributed in the range $0\le H<2M$.

@$@<If T is not in S, insert it and increase size@>==@{@-
H := trunc(alpha * (T-1));
while hash[H] > T do
  if H = 0 then H := H_max else H := H-1;
if hash[H] < T then { T is not present }
  begin size := size + 1;
    @<Insert T into the ordered hash table@>;
  end@}

@!------------------------------------------------------------

The heart of ordered hashing is the insertion process.  In general,
the new key @{T@} will be inserted in place of a previous key $T_1<T$,
which is then re-inserted in place of $T_2<T_1$, etc., until an empty
slot is discovered.

@$@<Insert T into the ordered hash table@>==@{@-
while hash[H] > 0 do
  begin TT := hash[H]; { we have 0 < TT < T }
  hash[H] := T; T := TT;
  repeat if H = 0 then H := H_max
    else H := H - 1;
  until hash[H] < T;
end;
hash[H] := T@}

@!------------------------------------------------------------
@a

@$@<Global variables@>+=@{@-
TT: integer; { a key that's being moved }
@}

@!------------------------------------------------------------
@a@<Sorting in linear time@>

The climax of this program is the fact that the entries in our ordered
hash table can easily be read out in increasing order.

Why is this true?  Well, we know that the final state of the table is
independent of the order in which the elements entered.  Furthermore
it's easy to understand what the table looks like when the entries are
inserted in decreasing order, because we have used a monotonic hash
function.  Therefore we know that the table must have an especially
simple form.

Suppose the nonzero entries are $T_1<\cdots<T_M$.  If $k$ of these
have `wrapped around' in the insertion process (i.e., if @{H@} passed
from 0 to @{H\_max@}, $k$ times), table position @{hash[0]@} will
either be zero (in which case $k$ must also be zero) or it will
contain $T_{k+1}$.  In the latter case, the entries
$T_{k+1}<\cdots<T_M$ and $T_1<\cdots<T_k$ will appear in order from
left to right.  Thus the output can be sorted with at most two passes
over the table!

@$@<print_it@>@M==@{@<print_ln@>@(hash[H] : 10@)@}

@$@<Print the elements of S in sorted order@>==@{@-
if hash[0] = 0 then { there was no wrap-around }
  begin for H := 1 to H_max do
    if hash[H] > 0 then @<print_it@>;
  end
else begin for H := 1 to H_max do
         { print the wrapped-around entries }
    if hash[H] > 0 then
      if hash[H] < hash[0] then @<print_it@>;
  for H := 0 to H_max do
    if hash[H] >= hash[0] then @<print_it@>;
  end@}