Commit | Line | Data |
---|---|---|
7f918cf1 CE |
1 | (* From the SML/NJ benchmark suite. *) |
2 | (* tree.sml | |
3 | * | |
4 | * COPYRIGHT (c) 1994 AT&T Bell Laboratories. | |
5 | * | |
6 | * Trees for the TSP program. | |
7 | *) | |
8 | ||
9 | structure Tree = | |
10 | struct | |
11 | ||
12 | datatype tree | |
13 | = NULL | |
14 | | ND of { | |
15 | left : tree, right : tree, | |
16 | x : real, y : real, | |
17 | sz : int, | |
18 | prev : tree ref, next : tree ref | |
19 | } | |
20 | ||
21 | fun mkNode (l, r, x, y, sz) = ND{ | |
22 | left = l, right = r, x = x, y = y, sz = sz, | |
23 | prev = ref NULL, next = ref NULL | |
24 | } | |
25 | ||
26 | fun printTree (outS, NULL) = () | |
27 | | printTree (outS, ND{x, y, left, right, ...}) = ( | |
28 | TextIO.output(outS, String.concat [ | |
29 | Real.toString x, " ", Real.toString y, "\n"]); | |
30 | printTree (outS, left); | |
31 | printTree (outS, right)) | |
32 | ||
33 | fun printList (outS, NULL) = () | |
34 | | printList (outS, start as ND{next, ...}) = let | |
35 | fun cycle (ND{next=next', ...}) = (next = next') | |
36 | | cycle _ = false | |
37 | fun prt (NULL) = () | |
38 | | prt (t as ND{x, y, next, ...}) = ( | |
39 | TextIO.output(outS, String.concat [ | |
40 | Real.toString x, " ", Real.toString y, "\n" | |
41 | ]); | |
42 | if (cycle (!next)) | |
43 | then () | |
44 | else prt (!next)) | |
45 | in | |
46 | prt start | |
47 | end | |
48 | ||
49 | end; | |
50 | ||
51 | (* tsp.sml | |
52 | * | |
53 | * COPYRIGHT (c) 1994 AT&T Bell Laboratories. | |
54 | *) | |
55 | ||
56 | structure TSP : sig | |
57 | ||
58 | val tsp : (Tree.tree * int) -> Tree.tree | |
59 | ||
60 | end = struct | |
61 | ||
62 | structure T = Tree | |
63 | ||
64 | fun setPrev (T.ND{prev, ...}, x) = prev := x | |
65 | fun setNext (T.ND{next, ...}, x) = next := x | |
66 | fun link (a as T.ND{next, ...}, b as T.ND{prev, ...}) = ( | |
67 | next := b; prev := a) | |
68 | ||
69 | fun sameNd (T.ND{next, ...}, T.ND{next=next', ...}) = (next = next') | |
70 | | sameNd (T.NULL, T.NULL) = true | |
71 | | sameNd _ = false | |
72 | ||
73 | (* Find Euclidean distance from a to b *) | |
74 | fun distance (T.ND{x=ax, y=ay, ...}, T.ND{x=bx, y=by, ...}) = | |
75 | Math.sqrt(((ax-bx)*(ax-bx)+(ay-by)*(ay-by))) | |
76 | | distance _ = raise Fail "distance" | |
77 | ||
78 | (* sling tree nodes into a list -- requires root to be tail of list, and | |
79 | * only fills in next field, not prev. | |
80 | *) | |
81 | fun makeList T.NULL = T.NULL | |
82 | | makeList (t as T.ND{left, right, next = t_next, ...}) = let | |
83 | val retVal = (case (makeList left, makeList right) | |
84 | of (T.NULL, T.NULL) => t | |
85 | | (l as T.ND{...}, T.NULL) => (setNext(left, t); l) | |
86 | | (T.NULL, r as T.ND{...}) => (setNext(right, t); r) | |
87 | | (l as T.ND{...}, r as T.ND{...}) => ( | |
88 | setNext(right, t); setNext(left, r); l) | |
89 | (* end case *)) | |
90 | in | |
91 | t_next := T.NULL; | |
92 | retVal | |
93 | end | |
94 | ||
95 | (* reverse orientation of list *) | |
96 | fun reverse T.NULL = () | |
97 | | reverse (t as T.ND{next, prev, ...}) = let | |
98 | fun rev (_, T.NULL) = () | |
99 | | rev (back, tmp as T.ND{prev, next, ...}) = let | |
100 | val tmp' = !next | |
101 | in | |
102 | next := back; setPrev(back, tmp); | |
103 | rev (tmp, tmp') | |
104 | end | |
105 | in | |
106 | setNext (!prev, T.NULL); | |
107 | prev := T.NULL; | |
108 | rev (t, !next) | |
109 | end | |
110 | ||
111 | (* Use closest-point heuristic from Cormen Leiserson and Rivest *) | |
112 | fun conquer (T.NULL) = T.NULL | |
113 | | conquer t = let | |
114 | val (cycle as T.ND{next=cycle_next, prev=cycle_prev, ...}) = makeList t | |
115 | fun loop (T.NULL) = () | |
116 | | loop (t as T.ND{next=ref doNext, prev, ...}) = | |
117 | let | |
118 | fun findMinDist (min, minDist, tmp as T.ND{next, ...}) = | |
119 | if (sameNd(cycle, tmp)) | |
120 | then min | |
121 | else let | |
122 | val test = distance(t, tmp) | |
123 | in | |
124 | if (test < minDist) | |
125 | then findMinDist (tmp, test, !next) | |
126 | else findMinDist (min, minDist, !next) | |
127 | end | |
128 | val (min as T.ND{next=ref min_next, prev=ref min_prev, ...}) = | |
129 | findMinDist (cycle, distance(t, cycle), !cycle_next) | |
130 | val minToNext = distance(min, min_next) | |
131 | val minToPrev = distance(min, min_prev) | |
132 | val tToNext = distance(t, min_next) | |
133 | val tToPrev = distance(t, min_prev) | |
134 | in | |
135 | if ((tToPrev - minToPrev) < (tToNext - minToNext)) | |
136 | then ( (* insert between min and min_prev *) | |
137 | link (min_prev, t); | |
138 | link (t, min)) | |
139 | else ( | |
140 | link (min, t); | |
141 | link (t, min_next)); | |
142 | loop doNext | |
143 | end | |
144 | val t' = !cycle_next | |
145 | in | |
146 | (* Create initial cycle *) | |
147 | cycle_next := cycle; cycle_prev := cycle; | |
148 | loop t'; | |
149 | cycle | |
150 | end | |
151 | ||
152 | (* Merge two cycles as per Karp *) | |
153 | fun merge (a as T.ND{next, ...}, b, t) = let | |
154 | fun locateCycle (start as T.ND{next, ...}) = let | |
155 | fun findMin (min, minDist, tmp as T.ND{next, ...}) = | |
156 | if (sameNd(start, tmp)) | |
157 | then (min, minDist) | |
158 | else let val test = distance(t, tmp) | |
159 | in | |
160 | if (test < minDist) | |
161 | then findMin (tmp, test, !next) | |
162 | else findMin (min, minDist, !next) | |
163 | end | |
164 | val (min as T.ND{next=ref next', prev=ref prev', ...}, minDist) = | |
165 | findMin (start, distance(t, start), !next) | |
166 | val minToNext = distance(min, next') | |
167 | val minToPrev = distance(min, prev') | |
168 | val tToNext = distance(t, next') | |
169 | val tToPrev = distance(t, prev') | |
170 | in | |
171 | if ((tToPrev - minToPrev) < (tToNext - minToNext)) | |
172 | (* would insert between min and prev *) | |
173 | then (prev', tToPrev, min, minDist) | |
174 | (* would insert between min and next *) | |
175 | else (min, minDist, next', tToNext) | |
176 | end | |
177 | (* Compute location for first cycle *) | |
178 | val (p1, tToP1, n1, tToN1) = locateCycle a | |
179 | (* compute location for second cycle *) | |
180 | val (p2, tToP2, n2, tToN2) = locateCycle b | |
181 | (* Now we have 4 choices to complete: | |
182 | * 1:t,p1 t,p2 n1,n2 | |
183 | * 2:t,p1 t,n2 n1,p2 | |
184 | * 3:t,n1 t,p2 p1,n2 | |
185 | * 4:t,n1 t,n2 p1,p2 | |
186 | *) | |
187 | val n1ToN2 = distance(n1, n2) | |
188 | val n1ToP2 = distance(n1, p2) | |
189 | val p1ToN2 = distance(p1, n2) | |
190 | val p1ToP2 = distance(p1, p2) | |
191 | fun choose (testChoice, test, choice, minDist) = | |
192 | if (test < minDist) then (testChoice, test) else (choice, minDist) | |
193 | val (choice, minDist) = (1, tToP1+tToP2+n1ToN2) | |
194 | val (choice, minDist) = choose(2, tToP1+tToN2+n1ToP2, choice, minDist) | |
195 | val (choice, minDist) = choose(3, tToN1+tToP2+p1ToN2, choice, minDist) | |
196 | val (choice, minDist) = choose(4, tToN1+tToN2+p1ToP2, choice, minDist) | |
197 | in | |
198 | case choice | |
199 | of 1 => ( (* 1:p1,t t,p2 n2,n1 -- reverse 2! *) | |
200 | reverse n2; | |
201 | link (p1, t); | |
202 | link (t, p2); | |
203 | link (n2, n1)) | |
204 | | 2 => ( (* 2:p1,t t,n2 p2,n1 -- OK *) | |
205 | link (p1, t); | |
206 | link (t, n2); | |
207 | link (p2, n1)) | |
208 | | 3 => ( (* 3:p2,t t,n1 p1,n2 -- OK *) | |
209 | link (p2, t); | |
210 | link (t, n1); | |
211 | link (p1, n2)) | |
212 | | 4 => ( (* 4:n1,t t,n2 p2,p1 -- reverse 1! *) | |
213 | reverse n1; | |
214 | link (n1, t); | |
215 | link (t, n2); | |
216 | link (p2, p1)) | |
217 | (* end case *); | |
218 | t | |
219 | end (* merge *) | |
220 | ||
221 | (* Compute TSP for the tree t -- use conquer for problems <= sz * *) | |
222 | fun tsp (t as T.ND{left, right, sz=sz', ...}, sz) = | |
223 | if (sz' <= sz) | |
224 | then conquer t | |
225 | else merge (tsp(left, sz), tsp(right, sz), t) | |
226 | | tsp (T.NULL, _) = T.NULL | |
227 | ||
228 | end; | |
229 | ||
230 | (* rand-sig.sml | |
231 | * | |
232 | * COPYRIGHT (c) 1993 by AT&T Bell Laboratories. See COPYRIGHT file for details. | |
233 | * COPYRIGHT (c) 1998 by AT&T Laboratories. | |
234 | * | |
235 | * Signature for a simple random number generator. | |
236 | * | |
237 | *) | |
238 | ||
239 | signature RAND = | |
240 | sig | |
241 | ||
242 | type rand = Word.word | |
243 | ||
244 | val randMin : rand | |
245 | val randMax : rand | |
246 | ||
247 | val random : rand -> rand | |
248 | (* Given seed, return value randMin <= v <= randMax | |
249 | * Iteratively using the value returned by random as the | |
250 | * next seed to random will produce a sequence of pseudo-random | |
251 | * numbers. | |
252 | *) | |
253 | ||
254 | val mkRandom : rand -> unit -> rand | |
255 | (* Given seed, return function generating a sequence of | |
256 | * random numbers randMin <= v <= randMax | |
257 | *) | |
258 | ||
259 | val norm : rand -> real | |
260 | (* Map values in the range [randMin,randMax] to (0.0,1.0) *) | |
261 | ||
262 | val range : (int * int) -> rand -> int | |
263 | (* Map v, randMin <= v <= randMax, to integer range [i,j] | |
264 | * Exception - | |
265 | * Fail if j < i | |
266 | *) | |
267 | ||
268 | end (* RAND *) | |
269 | ||
270 | (* rand.sml | |
271 | * | |
272 | * COPYRIGHT (c) 1991 by AT&T Bell Laboratories. See COPYRIGHT file for details | |
273 | * COPYRIGHT (c) 1998 by AT&T Laboratories. See COPYRIGHT file for details | |
274 | * | |
275 | * Random number generator taken from Paulson, pp. 170-171. | |
276 | * Recommended by Stephen K. Park and Keith W. Miller, | |
277 | * Random number generators: good ones are hard to find, | |
278 | * CACM 31 (1988), 1192-1201 | |
279 | * Updated to include the new preferred multiplier of 48271 | |
280 | * CACM 36 (1993), 105-110 | |
281 | * Updated to use on Word. | |
282 | * | |
283 | * Note: The Random structure provides a better generator. | |
284 | *) | |
285 | ||
286 | structure Rand : RAND = | |
287 | struct | |
288 | ||
289 | type rand = Word.word | |
290 | type rand' = Int.int (* internal representation *) | |
291 | ||
292 | val a : rand' = 48271 | |
293 | val m : rand' = valOf Int.maxInt (* 2^31 - 1 *) | |
294 | val m_1 = m - 1 | |
295 | val q = m div a | |
296 | val r = m mod a | |
297 | ||
298 | val extToInt = Word.toInt | |
299 | val intToExt = Word.fromInt | |
300 | ||
301 | val randMin : rand = 0w1 | |
302 | val randMax : rand = intToExt m_1 | |
303 | ||
304 | fun chk 0w0 = 1 | |
305 | | chk 0wx7fffffff = m_1 | |
306 | | chk seed = extToInt seed | |
307 | ||
308 | fun random' seed = let | |
309 | val hi = seed div q | |
310 | val lo = seed mod q | |
311 | val test = a * lo - r * hi | |
312 | in | |
313 | if test > 0 then test else test + m | |
314 | end | |
315 | ||
316 | val random = intToExt o random' o chk | |
317 | ||
318 | fun mkRandom seed = let | |
319 | val seed = ref (chk seed) | |
320 | in | |
321 | fn () => (seed := random' (!seed); intToExt (!seed)) | |
322 | end | |
323 | ||
324 | val real_m = Real.fromInt m | |
325 | fun norm s = (Real.fromInt (Word.toInt s)) / real_m | |
326 | ||
327 | fun range (i,j) = | |
328 | if j < i | |
329 | then raise Fail "Random.range: hi < lo" | |
330 | else if j = i then fn _ => i | |
331 | else let | |
332 | val R = Int.fromInt j - Int.fromInt i | |
333 | val cvt = Word.toIntX o Word.fromInt | |
334 | in | |
335 | if R = m then Word.toIntX | |
336 | else fn s => i + cvt ((extToInt s) mod (R+1)) | |
337 | end | |
338 | ||
339 | end (* Rand *) | |
340 | ||
341 | (* build.sml | |
342 | * | |
343 | * COPYRIGHT (c) 1994 AT&T Bell Laboratories. | |
344 | * | |
345 | * Build a two-dimensional tree for TSP. | |
346 | *) | |
347 | ||
348 | structure BuildTree : sig | |
349 | ||
350 | datatype axis = X_AXIS | Y_AXIS | |
351 | ||
352 | val buildTree : { | |
353 | n : int, dir : axis, | |
354 | min_x : real, min_y : real, max_x : real, max_y : real | |
355 | } -> Tree.tree | |
356 | ||
357 | end = struct | |
358 | ||
359 | structure T = Tree | |
360 | ||
361 | val m_e = 2.7182818284590452354 | |
362 | val m_e2 = 7.3890560989306502274 | |
363 | val m_e3 = 20.08553692318766774179 | |
364 | val m_e6 = 403.42879349273512264299 | |
365 | val m_e12 = 162754.79141900392083592475 | |
366 | ||
367 | datatype axis = X_AXIS | Y_AXIS | |
368 | ||
369 | (* builds a 2D tree of n nodes in specified range with dir as primary axis *) | |
370 | fun buildTree arg = let | |
371 | val rand = Rand.mkRandom 0w314 | |
372 | fun drand48 () = Rand.norm (rand ()) | |
373 | fun median {min, max, n} = let | |
374 | val t = drand48(); (* in [0.0..1.0) *) | |
375 | val retval = if (t > 0.5) | |
376 | then Math.ln(1.0-(2.0*(m_e12-1.0)*(t-0.5)/m_e12))/12.0 | |
377 | else ~(Math.ln(1.0-(2.0*(m_e12-1.0)*t/m_e12))/12.0) | |
378 | in | |
379 | min + ((retval + 1.0) * (max - min)/2.0) | |
380 | end | |
381 | fun uniform {min, max} = min + (drand48() * (max - min)) | |
382 | fun build {n = 0, ...} = T.NULL | |
383 | | build {n, dir=X_AXIS, min_x, min_y, max_x, max_y} = let | |
384 | val med = median{min=min_y, max=max_y, n=n} | |
385 | fun mkTree (min, max) = build{ | |
386 | n=n div 2, dir=Y_AXIS, min_x=min_x, max_x=max_x, | |
387 | min_y=min, max_y=max | |
388 | } | |
389 | in | |
390 | T.mkNode( | |
391 | mkTree(min_y, med), mkTree(med, max_y), | |
392 | uniform{min=min_x, max=max_x}, med, n) | |
393 | end | |
394 | | build {n, dir=Y_AXIS, min_x, min_y, max_x, max_y} = let | |
395 | val med = median{min=min_x, max=max_x, n=n} | |
396 | fun mkTree (min, max) = build{ | |
397 | n=n div 2, dir=X_AXIS, min_x=min, max_x=max, | |
398 | min_y=min_y, max_y=max_y | |
399 | } | |
400 | in | |
401 | T.mkNode( | |
402 | mkTree(min_x, med), mkTree(med, max_x), | |
403 | med, uniform{min=min_y, max=max_y}, n) | |
404 | end | |
405 | in | |
406 | build arg | |
407 | end | |
408 | ||
409 | end; (* Build *) | |
410 | ||
411 | signature BMARK = | |
412 | sig | |
413 | val doit : int -> unit | |
414 | val testit : TextIO.outstream -> unit | |
415 | end; | |
416 | (* main.sml | |
417 | * | |
418 | * COPYRIGHT (c) 1994 AT&T Bell Laboratories. | |
419 | * | |
420 | *) | |
421 | ||
422 | structure Main : sig | |
423 | ||
424 | include BMARK | |
425 | ||
426 | val dumpPS : TextIO.outstream -> unit | |
427 | ||
428 | end = struct | |
429 | ||
430 | val name = "TSP" | |
431 | ||
432 | val problemSz = ref 32767 | |
433 | val divideSz = ref 150 | |
434 | ||
435 | fun printLength (outS, Tree.NULL) = print "(* 0 points *)\n" | |
436 | | printLength (outS, start as Tree.ND{next, x, y, ...}) = let | |
437 | fun cycle (Tree.ND{next=next', ...}) = (next = next') | |
438 | | cycle _ = false | |
439 | fun distance (ax, ay, bx, by) = let | |
440 | val dx = ax-bx and dy = ay-by | |
441 | in | |
442 | Math.sqrt (dx*dx + dy*dy) | |
443 | end | |
444 | fun length (Tree.NULL, px, py, n, len) = (n, len+distance(px, py, x, y)) | |
445 | | length (t as Tree.ND{x, y, next, ...}, px, py, n, len) = | |
446 | if (cycle t) | |
447 | then (n, len+distance(px, py, x, y)) | |
448 | else length(!next, x, y, n+1, len+distance(px, py, x, y)) | |
449 | in | |
450 | if (cycle(!next)) | |
451 | then TextIO.output (outS, "(* 1 point *)\n") | |
452 | else let | |
453 | val (n, len) = length(!next, x, y, 1, 0.0) | |
454 | in | |
455 | TextIO.output (outS, concat[ | |
456 | "(* ", Int.toString n, "points, cycle length = ", | |
457 | Real.toString len, " *)\n" | |
458 | ]) | |
459 | end | |
460 | end | |
461 | ||
462 | fun mkTree n = BuildTree.buildTree { | |
463 | n=n, dir=BuildTree.X_AXIS, | |
464 | min_x=0.0, max_x=1.0, | |
465 | min_y=0.0, max_y=1.0 | |
466 | } | |
467 | ||
468 | fun doit' n = TSP.tsp (mkTree n, !divideSz) | |
469 | ||
470 | fun dumpPS outS = ( | |
471 | TextIO.output (outS, "newgraph\n"); | |
472 | TextIO.output (outS, "newcurve pts\n"); | |
473 | Tree.printList (outS, doit' (!problemSz)); | |
474 | TextIO.output (outS, "linetype solid\n")) | |
475 | ||
476 | fun testit strm = printLength (strm, doit' (!problemSz)) | |
477 | ||
478 | val _ = problemSz := 2097151 | |
479 | fun doit () = doit' (!problemSz) | |
480 | ||
481 | val doit = | |
482 | fn n => | |
483 | let | |
484 | fun loop n = | |
485 | if n = 0 | |
486 | then () | |
487 | else (doit(); | |
488 | loop(n-1)) | |
489 | in loop n | |
490 | end | |
491 | ||
492 | end |