1 /***************************************************************************/
2 /* */
3 /* FUNCTIONS FOR COMPUTING EDGE LENGTHS FOR GEOMETRIC PROBLEMS */
4 /* */
5 /* TSP CODE */
6 /* */
7 /* */
8 /* Written by: Applegate, Bixby, Chvatal, and Cook */
9 /* Date: Summer 1994 */
10 /* Modified - March 2, 1995 */
11 /* - October 5, 1995 (Bico) */
12 /* */
13 /* */
14 /* EXPORTED FUNCTIONS: */
15 /* int */
16 /* CCutil_init_dat_edgelen (CCdatagroup *dat) */
17 /* int */
18 /* (*CCutil_dat_edgelen) (int i, int j, CCdatagroup *dat) */
19 /* void */
20 /* CCutil_dsjrand_init (int maxdist, int seed) */
21 /* void */
22 /* CCutil_freedata (int ncount, CCdatagroup *dat) */
23 /* */
24 /* NOTES: */
25 /* Supported norms (with defs in edgelen.h) are: */
26 /* MAXNORM - the L-infinity norm */
27 /* EUCLIDEAN_CEIL - the norm for the plaXXXX problems */
28 /* EUCLIDEAN - rounded L-2 norm */
29 /* EUCLIDEAN_3d - rounded L-2 norm in 3 space */
30 /* IBM - a norm for drilling problems in Bonn */
31 /* GEOGRAPHIC - distances on a sphere (Groetshel and Holland) */
32 /* ATT - pseudo-Euclidean norm for att532 */
33 /* MATRIXNORM - complete graph (lower + diagonal matrix) */
34 /* DSJRAND - random edgelengths */
35 /* CRYSTAL - Bland-Shallcross xray norm */
36 /* - The coordinates generated for CRYSTAL problems */
37 /* (in getdata.c) have been diveded by the motor speeds */
38 /* (this makes the edgelen function faster) and scaled by */
39 /* CRYSTAL_SCALE (currently 10000) and rounded to the */
40 /* nearest integer (this lets the edgelen function produce */
41 /* integer lengths without further rounding). The result is */
42 /* a closer approximation to the Bland - Shallcrosss */
43 /* floating point length function than that given in */
44 /* TSPLIB_1.2. */
45 /* */
46 /* */
47 /***************************************************************************/
48
49 #include "machdefs.h"
50 #include "util.h"
51
52 #ifdef CC_PROTOTYPE_ANSI
53
54 static double
55 dtrunc (double);
56
57 #else
58
59 static double
60 dtrunc ();
61
62 #endif
63
64 #ifndef M_PI
65 #define M_PI 3.14159265358979323846264
66 #endif
67
68 #ifdef CC_PROTOTYPE_ANSI
69 int (*CCutil_dat_edgelen) (int, int, CCdatagroup *);
70 #else
71 int (*CCutil_dat_edgelen) ();
72 #endif
73
74 #ifdef CC_PROTOTYPE_ANSI
CCutil_init_dat_edgelen(CCdatagroup * dat)75 int CCutil_init_dat_edgelen (CCdatagroup *dat)
76 #else
77 int CCutil_init_dat_edgelen (dat)
78 CCdatagroup *dat;
79 #endif
80 {
81 switch (dat->norm) {
82 case CC_EUCLIDEAN_CEIL:
83 CCutil_dat_edgelen = CCutil_euclid_ceiling_edgelen;
84 break;
85 case CC_EUCLIDEAN:
86 CCutil_dat_edgelen = CCutil_euclid_edgelen;
87 break;
88 case CC_MAXNORM:
89 CCutil_dat_edgelen = CCutil_max_edgelen;
90 break;
91 case CC_EUCLIDEAN_3D:
92 CCutil_dat_edgelen = CCutil_euclid3d_edgelen;
93 break;
94 case CC_IBM:
95 CCutil_dat_edgelen = CCutil_ibm_edgelen;
96 break;
97 case CC_GEOGRAPHIC:
98 CCutil_dat_edgelen = CCutil_geographic_edgelen;
99 break;
100 case CC_ATT:
101 CCutil_dat_edgelen = CCutil_att_edgelen;
102 break;
103 case CC_MATRIXNORM:
104 CCutil_dat_edgelen = CCutil_matrix_edgelen;
105 break;
106 case CC_DSJRANDNORM:
107 CCutil_dat_edgelen = CCutil_dsjrand_edgelen;
108 break;
109 case CC_CRYSTAL:
110 CCutil_dat_edgelen = CCutil_crystal_edgelen;
111 break;
112 default:
113 fprintf (stderr, "ERROR: Unknown NORM %d.\n", dat->norm);
114 return 1;
115 }
116 return 0;
117 }
118
119 /* Several variables that would normally be called y1 and y2 are called
120 yy1 and yyy2 to avoid conflict with the bessel functions */
121
122 #ifdef CC_PROTOTYPE_ANSI
CCutil_max_edgelen(int i,int j,CCdatagroup * dat)123 int CCutil_max_edgelen (int i, int j, CCdatagroup *dat)
124 #else
125 int CCutil_max_edgelen (i, j, dat)
126 int i, j;
127 CCdatagroup *dat;
128 #endif
129 {
130 double t1 = dat->x[i] - dat->x[j], t2 = dat->y[i] - dat->y[j];
131
132 if (t1 < 0)
133 t1 *= -1;
134 if (t2 < 0)
135 t2 *= -1;
136
137 return (int) (t1 < t2 ? t2 : t1);
138 }
139
140 #ifdef CC_PROTOTYPE_ANSI
CCutil_euclid_edgelen(int i,int j,CCdatagroup * dat)141 int CCutil_euclid_edgelen (int i, int j, CCdatagroup *dat)
142 #else
143 int CCutil_euclid_edgelen (i, j, dat)
144 int i, j;
145 CCdatagroup *dat;
146 #endif
147 {
148 double t1 = dat->x[i] - dat->x[j], t2 = dat->y[i] - dat->y[j];
149 int temp;
150
151 temp = (int) (sqrt (t1 * t1 + t2 * t2) + 0.5);
152 return temp;
153 }
154
155 #ifdef CC_PROTOTYPE_ANSI
CCutil_euclid3d_edgelen(int i,int j,CCdatagroup * dat)156 int CCutil_euclid3d_edgelen (int i, int j, CCdatagroup *dat)
157 #else
158 int CCutil_euclid3d_edgelen (i, j, dat)
159 int i, j;
160 CCdatagroup *dat;
161 #endif
162 {
163 double t1 = dat->x[i] - dat->x[j], t2 = dat->y[i] - dat->y[j];
164 double t3 = dat->z[i] - dat->z[j];
165 int temp;
166
167 temp = (int) (sqrt (t1 * t1 + t2 * t2 + t3 * t3) + 0.5);
168 return temp;
169 }
170
171 #ifdef CC_PROTOTYPE_ANSI
CCutil_ibm_edgelen(int i,int j,CCdatagroup * dat)172 int CCutil_ibm_edgelen (int i, int j, CCdatagroup *dat)
173 #else
174 int CCutil_ibm_edgelen (i, j, dat)
175 int i, j;
176 CCdatagroup *dat;
177 #endif
178 {
179 double dw = dat->x[i] - dat->x[j], dw1 = dat->y[i] - dat->y[j];
180 static double ibm_xmult[7] = {1062.5,
181 300.0,
182 300.0,
183 250.0,
184 300.0,
185 1000.0,
186 154.6};
187 static double ibm_xadd[7] = {155.0 - 0.01 * 1062.5,
188 197.5 - 0.05 * 300.0,
189 212.5 - 0.10 * 300.0,
190 227.5 - 0.15 * 250.0,
191 240.5 - 0.20 * 300.0,
192 255.0 - 0.25 * 1000.0,
193 305.0 - 0.30 * 154.6};
194 static double ibm_ymult[7] = {1062.5,
195 450.0,
196 350.0,
197 250.0,
198 300.0,
199 900.0,
200 157.7};
201 static double ibm_yadd[7] = {150.0 - 0.01 * 1062.5,
202 192.5 - 0.05 * 450.0,
203 215.0 - 0.10 * 350.0,
204 232.5 - 0.15 * 250.0,
205 245.5 - 0.20 * 300.0,
206 250.0 - 0.25 * 900.0,
207 295.0 - 0.30 * 157.7};
208
209 if (dw < 0.0)
210 dw = -dw;
211 dw /= 25400.0;
212 if (dw <= 0.01) {
213 dw *= 15500.0;
214 } else if (dw >= 0.30) {
215 dw = dw * 154.6 + (305.0 - 0.3 * 154.6);
216 } else {
217 dw = dw * ibm_xmult[(int) (dw / 0.05)] +
218 ibm_xadd[(int) (dw / 0.05)];
219 }
220 if (dw1 < 0.0)
221 dw1 = -dw1;
222 dw1 /= 25400.0;
223 if (dw1 <= 0.01) {
224 dw1 *= 15000.0;
225 } else if (dw1 >= 0.30) {
226 dw1 = dw1 * 157.7 + (295.0 - 0.3 * 157.7);
227 } else {
228 dw1 = dw1 * ibm_ymult[(int) (dw1 / 0.05)] +
229 ibm_yadd[(int) (dw1 / 0.05)];
230 }
231 if (dw < dw1)
232 dw = dw1;
233 return (int) dw;
234 }
235
236 #ifdef CC_PROTOTYPE_ANSI
CCutil_euclid_ceiling_edgelen(int i,int j,CCdatagroup * dat)237 int CCutil_euclid_ceiling_edgelen (int i, int j, CCdatagroup *dat)
238 #else
239 int CCutil_euclid_ceiling_edgelen (i, j, dat)
240 int i, j;
241 CCdatagroup *dat;
242 #endif
243 {
244 double t1 = dat->x[i] - dat->x[j], t2 = dat->y[i] - dat->y[j];
245 /*
246 int rd;
247 double max;
248
249 max = sqrt (t1 * t1 + t2 * t2);
250 rd = (int) max;
251 return (((max - rd) > .000000001) ? rd + 1 : rd);
252 */
253 return (int) (ceil (sqrt (t1 * t1 + t2 * t2)));
254 }
255
256 #ifdef CC_PROTOTYPE_ANSI
CCutil_geographic_edgelen(int i,int j,CCdatagroup * dat)257 int CCutil_geographic_edgelen (int i, int j, CCdatagroup *dat)
258 #else
259 int CCutil_geographic_edgelen (i, j, dat)
260 int i, j;
261 CCdatagroup *dat;
262 #endif
263 {
264 double deg, min;
265 double lati, latj, longi, longj;
266 double q1, q2, q3;
267 int dd;
268 double x1 = dat->x[i], x2 = dat->x[j], yy1 = dat->y[i], yy2 = dat->y[j];
269
270 deg = dtrunc (x1);
271 min = x1 - deg;
272 lati = M_PI * (deg + 5.0 * min / 3.0) / 180.0;
273 deg = dtrunc (x2);
274 min = x2 - deg;
275 latj = M_PI * (deg + 5.0 * min / 3.0) / 180.0;
276
277 deg = dtrunc (yy1);
278 min = yy1 - deg;
279 longi = M_PI * (deg + 5.0 * min / 3.0) / 180.0;
280 deg = dtrunc (yy2);
281 min = yy2 - deg;
282 longj = M_PI * (deg + 5.0 * min / 3.0) / 180.0;
283
284 q1 = cos (longi - longj);
285 q2 = cos (lati - latj);
286 q3 = cos (lati + latj);
287 dd = (int) (6378.388 * acos (0.5 * ((1.0 + q1) * q2 - (1.0 - q1) * q3))
288 + 1.0);
289 return dd;
290 }
291
292 #ifdef CC_PROTOTYPE_ANSI
CCutil_att_edgelen(int i,int j,CCdatagroup * dat)293 int CCutil_att_edgelen (int i, int j, CCdatagroup *dat)
294 #else
295 int CCutil_att_edgelen (i, j, dat)
296 int i, j;
297 CCdatagroup *dat;
298 #endif
299 {
300 double xd = dat->x[i] - dat->x[j];
301 double yd = dat->y[i] - dat->y[j];
302 double rij = sqrt ((xd * xd + yd * yd) / 10.0);
303 double tij = dtrunc (rij);
304 int dij;
305
306 if (tij < rij)
307 dij = (int) tij + 1;
308 else
309 dij = (int) tij;
310 return dij;
311 }
312
313 #ifdef CC_PROTOTYPE_ANSI
dtrunc(double x)314 static double dtrunc (double x)
315 #else
316 static double dtrunc (x)
317 double x;
318 #endif
319 {
320 int k;
321
322 k = (int) x;
323 x = (double) k;
324 return x;
325 }
326
327 static int dsjrand_param = 1;
328 static double dsjrand_factor = 1.0;
329
330 #ifdef CC_PROTOTYPE_ANSI
CCutil_dsjrand_init(int maxdist,int seed)331 void CCutil_dsjrand_init (int maxdist, int seed)
332 #else
333 void CCutil_dsjrand_init (maxdist, seed)
334 int maxdist, seed;
335 #endif
336 {
337 dsjrand_factor = maxdist/2147483648.0;
338 dsjrand_param = 104*seed+1;
339 }
340
341 #ifdef CC_PROTOTYPE_ANSI
CCutil_dsjrand_edgelen(int i,int j,CCdatagroup * dat)342 int CCutil_dsjrand_edgelen (int i, int j, CCdatagroup *dat)
343 #else
344 int CCutil_dsjrand_edgelen (i, j, dat)
345 int i, j;
346 CCdatagroup *dat;
347 #endif
348 {
349 int di = (int) dat->x[i];
350 int dj = (int) dat->x[j];
351 int x, y, z;
352
353 x = di&dj;
354 y = di|dj;
355 z = dsjrand_param;
356
357 x *= z;
358 y *= x;
359 z *= y;
360
361 z ^= dsjrand_param;
362
363 x *= z;
364 y *= x;
365 z *= y;
366
367 x = ((di+dj)^z)&0x7fffffff;
368 return (int)(x*dsjrand_factor);
369 }
370
371 #define CRYSTAL_SCALE 10000
372
373 #define CRYSTAL_FLIP_TOL ((180 * CRYSTAL_SCALE * 4) / 5)
374 #define CRYSTAL_NEEDS_FLIP(x) ((x) > (CRYSTAL_FLIP_TOL))
375 #define CRYSTAL_FLIP(x) ((2 * (CRYSTAL_FLIP_TOL)) - (x))
376
377 #ifdef CC_PROTOTYPE_ANSI
CCutil_crystal_edgelen(int i,int j,CCdatagroup * dat)378 int CCutil_crystal_edgelen (int i, int j, CCdatagroup *dat)
379 #else
380 int CCutil_crystal_edgelen (i, j, dat)
381 int i, j;
382 CCdatagroup *dat;
383 #endif
384 {
385 double w, w1;
386
387 w = dat->x[i] - dat->x[j];
388 if (w < 0)
389 w = -w;
390 w1 = dat->y[i] - dat->y[j];
391 if (w1 < 0)
392 w1 = -w1;
393 if (CRYSTAL_NEEDS_FLIP (w1))
394 w1 = CRYSTAL_FLIP (w1);
395 if (w < w1)
396 w = w1;
397 w1 = dat->z[i] - dat->z[j];
398 if (w1 < 0)
399 w1 = -w1;
400 if (w < w1)
401 w = w1;
402
403 return (int) w;
404 }
405
406 #ifdef CC_PROTOTYPE_ANSI
CCutil_matrix_edgelen(int i,int j,CCdatagroup * dat)407 int CCutil_matrix_edgelen (int i, int j, CCdatagroup *dat)
408 #else
409 int CCutil_matrix_edgelen (i, j, dat)
410 int i, j;
411 CCdatagroup *dat;
412 #endif
413 {
414 if (i > j)
415 return (dat->adj[i])[j];
416 else
417 return (dat->adj[j])[i];
418 }
419
420 #ifdef CC_PROTOTYPE_ANSI
CCutil_freedatagroup(int ncount,CCdatagroup * dat)421 void CCutil_freedatagroup (int ncount, CCdatagroup *dat)
422 #else
423 void CCutil_freedatagroup (ncount, dat)
424 int ncount;
425 CCdatagroup *dat;
426 #endif
427 {
428 int i;
429
430 CC_IFFREE (dat->x, double);
431 CC_IFFREE (dat->y, double);
432 CC_IFFREE (dat->z, double);
433 if (dat->adj) {
434 for (i = 0; i < ncount; i++) {
435 CC_IFFREE (dat->adj[i], int);
436 }
437 CC_FREE (dat->adj, int *);
438 }
439 }
440