efst.c

生成直角Steiner树的程序包

			" [-v N]"
			" <points-file\n",
			me);
	pp = &arg_doc [0];
	while ((p = *pp++) NE NULL) {
		(void) fprintf (stderr, "%s\n", p);
	}
	exit (1);
}

/*
 * Use the heapsort algorithm to sort the given terminals in increasing
 * order by the following keys:
 *
 *	1.	X coordinate
 *	2.	Y coordinate
 *	3.	index (i.e., position within input data)
 *
 * Of course, we do not move the points, but rather permute an array
 * of indexes into the points.
 */

	static
	int *
heapsort_x (

struct pset *		pts		/* IN - the terminals to sort */
)
{
int			i, i1, i2, j, k, n;
struct point *		p1;
struct point *		p2;
int *			index;

	n = pts -> n;

	index = NEWA (n, int);
	for (i = 0; i < n; i++) {
		index [i] = i;
	}

	/* Construct the heap via sift-downs, in O(n) time. */
	for (k = n >> 1; k >= 0; k--) {
		j = k;
		for (;;) {
			i = (j << 1) + 1;
			if (i + 1 < n) {
				/* Increment i (to right subchild of j) */
				/* if the right subchild is greater. */
				i1 = index [i];
				i2 = index [i + 1];
				p1 = &(pts -> a [i1]);
				p2 = &(pts -> a [i2]);
				if ((p2 -> x > p1 -> x) OR
				    ((p2 -> x EQ p1 -> x) AND
				     ((p2 -> y > p1 -> y) OR
				      ((p2 -> y EQ p1 -> y) AND
				       (i2 > i1))))) {
					++i;
				}
			}
			if (i >= n) {
				/* Hit bottom of heap, sift-down is done. */
				break;
			}
			i1 = index [j];
			i2 = index [i];
			p1 = &(pts -> a [i1]);
			p2 = &(pts -> a [i2]);
			if ((p1 -> x > p2 -> x) OR
			    ((p1 -> x EQ p2 -> x) AND
			     ((p1 -> y > p2 -> y) OR
			      ((p1 -> y EQ p2 -> y) AND
			       (i1 > i2))))) {
				/* Greatest child is smaller.  Sift-	*/
				/* down is done. */
				break;
			}
			/* Sift down and continue. */
			index [j] = i2;
			index [i] = i1;
			j = i;
		}
	}

	/* Now do actual sorting.  Exchange first/last and sift down. */
	while (n > 1) {
		/* Largest is at index [0], swap with index [n-1],	*/
		/* thereby putting it into final position.		*/
		--n;
		i = index [0];
		index [0] = index [n];
		index [n] = i;

		/* Now restore the heap by sifting index [0] down. */
		j = 0;
		for (;;) {
			i = (j << 1) + 1;
			if (i + 1 < n) {
				/* Increment i (to right subchild of j) */
				/* if the right subchild is greater. */
				i1 = index [i];
				i2 = index [i + 1];
				p1 = &(pts -> a [i1]);
				p2 = &(pts -> a [i2]);
				if ((p2 -> x > p1 -> x) OR
				    ((p2 -> x EQ p1 -> x) AND
				     ((p2 -> y > p1 -> y) OR
				      ((p2 -> y EQ p1 -> y) AND
				       (i2 > i1))))) {
					++i;
				}
			}
			if (i >= n) {
				/* Hit bottom of heap, sift-down is done. */
				break;
			}
			i1 = index [j];
			i2 = index [i];
			p1 = &(pts -> a [i1]);
			p2 = &(pts -> a [i2]);
			if ((p1 -> x > p2 -> x) OR
			    ((p1 -> x EQ p2 -> x) AND
			     ((p1 -> y > p2 -> y) OR
			      ((p1 -> y EQ p2 -> y) AND
			       (i1 > i2))))) {
				/* Greatest child is smaller.  Sift-	*/
				/* down is done. */
				break;
			}
			/* Sift down and continue. */
			index [j] = i2;
			index [i] = i1;
			j = i;
		}
	}

	return (index);
}

/*
 * This routine finds all pairs of points whose coordinates are exactly
 * identical.  This is a degeneracy that can cause successive algorithms
 * much heartburn, so we take care of them immediately by keeping only
 * the earliest (lowest pnum) point, and marking the others as duplicates.
 * We generate a partition of such terminals into subsets -- if terminals
 * A and B reside in the same subset then they have identical coordinates.
 * We refer to such a subset as a Duplicate Terminal Group.
 *
 * For each terminal group we retain ONLY THE FIRST member -- the terminal
 * map bits are turned OFF for all other members of a terminal group,
 * effectively eliminating those terminals from the problem.
 */

	static
	int
generate_duplicate_terminal_groups (

struct pset *		pts,	/* IN - original terminal set */
int *			xorder, /* IN - terminal numbers sorted by X coord */
int ***			grps	/* OUT - duplicate terminal groups */
)
{
int			i;
int			j;
int			n;
int			n_grps;
struct point *		p0;
struct point *		p1;
struct point *		p2;
int *			ip;
int **			real_ptrs;
int *			real_terms;
int **			ptrs;
int *			terms;

	n = pts -> n;

	n_grps = 0;

	ptrs	= NEWA (n + 1, int *);
	terms	= NEWA (n, int);

	ip = &terms [0];
	for (i = 1; i < n; ) {
		p0 = &(pts -> a [xorder [i - 1]]);
		p1 = &(pts -> a [xorder [i]]);

		if ((p0 -> y NE p1 -> y) OR (p0 -> x NE p1 -> x)) {
			/* Not identical. */
			++i;
			continue;
		}

		/* Terminals xorder[i-1] and xorder[i] are identical. */

		for (j = i + 1; j < n; j++) {
			p2 = &(pts -> a [xorder [j]]);
			if (p0 -> y NE p2 -> y) break;
			if (p0 -> x NE p2 -> x) break;
		}
		/* j now points to the first non-equal terminal */

		/* Make a new duplicate terminal group... */
		ptrs [n_grps++] = ip;
		*ip++ = xorder [i - 1];
		while (i < j) {
			*ip++ = xorder [i++];
		}

		/* Skip whole group of coincident points and continue */
	}
	ptrs [n_grps] = ip;

	if (n_grps <= 0) {
		*grps = NULL;
	}
	else {
		/* Transfer to permanent memory of proper size... */
		real_ptrs = NEWA (n_grps + 1, int *);
		real_terms = NEWA (ip - terms, int);
		(void) memcpy ((char *) real_terms,
			       (char *) terms,
			       (ip - terms) * sizeof (int));
		for (i = 0; i <= n_grps; i++) {
			real_ptrs [i] = &real_terms [ptrs [i] - ptrs [0]];
		}
		*grps = real_ptrs;
	}

	free ((char *) terms);
	free ((char *) ptrs);

	return (n_grps);
}

/*
 * This routine removes all but the first terminal in each duplicate
 * terminal group.  We also prepare arrays for mapping forward from
 * old to new terminal numbers, and backward from new terminal numbers
 * to old.
 */

	static
	struct pset *
remove_duplicates (

struct pset *		pts,		/* IN - original point set */
int			ndg,		/* IN - number of duplicate groups */
int **			list,		/* IN - list of duplicate groups */
int **			fwd_map_ptr,	/* OUT - map to renumber old to new */
int **			rev_map_ptr	/* OUT - map to renumber new to old */
)
{
int		i;
int		j;
int		n;
int		kmasks;
int		numdel;
int		new_n;
int *		ip1;
int *		ip2;
int *		fwd;
int *		rev;
bitmap_t *	deleted;
struct pset *	newpts;

	n = pts -> n;

	kmasks = BMAP_ELTS (n);
	deleted = NEWA (kmasks, bitmap_t);
	for (i = 0; i < kmasks; i++) {
		deleted [i] = 0;
	}

	numdel = 0;
	for (i = 0; i < ndg; i++) {
		ip1 = list [i];
		ip2 = list [i + 1];

		/* Retain the first in this group, exclude all	*/
		/* of the others...				*/
		while (++ip1 < ip2) {
			if (BITON (deleted, *ip1)) {
				fatal ("remove_duplicates: Bug 1.");
			}
			++numdel;
			SETBIT (deleted, *ip1);
		}
	}

	new_n = n - numdel;

	fwd = NEWA (n, int);
	rev = NEWA (new_n, int);
	newpts = NEW_PSET (new_n);
	ZERO_PSET (newpts, new_n);
	newpts -> n = new_n;
	j = 0;
	for (i = 0; i < n; i++) {
		if (BITON (deleted, i)) {
			fwd [i] = -1;
		}
		else {
			newpts -> a [j].x		= pts -> a [i].x;
			newpts -> a [j].y		= pts -> a [i].y;
			newpts -> a [j].pnum		= j;
			rev [j] = i;
			fwd [i] = j;
			++j;
		}
	}

	free ((char *) deleted);

	*fwd_map_ptr = fwd;
	*rev_map_ptr = rev;

	return (newpts);
}


/* Solve a quadratic equation of the form
 * A/2 X^2 + B X + C/2 = 0
 *
 * We use a numerically robust formula from Press et al: Numerical
 * Recipies.
 * Only solves the equation if A <> 0.
 * Return FALSE if there are no real roots.
 */

	static
	bool
solve_quadratic (

double	 A,	 /* IN - twice the first coefficient */
double	 B,	 /* IN - second coefficient */
double	 C,	 /* IN - twice the third coefficient */
double * root1,	 /* OUT - first root */
double * root2	 /* OUT - second root */
)
{
	double Q, D;

	if (A EQ 0.0) return FALSE; /* this is not a quadratic equation */

	D = B*B - A*C;
	if (D < 0.0)  return FALSE; /* no real roots */

	D = sqrt(D);
	if (B >= 0.0)
		Q = - B - D;
	else
		Q = - B + D;

	if (Q NE 0.0) {
		*root1 = Q / A;
		*root2 = C / Q;
	}
	else {
		*root1 = 0.0;
		*root2 = -2.0 * B / A;
	}

	return TRUE;
}

/*
 * Upper bound heuristic (just calls the heuristic that was chosen)
 */

	static
	dist_t
upper_bound_heuristic (

struct pset *	pts,	/* IN - point set */
struct bsd  *	bsdp	/* IN - BSD info (can be NULL) */
)
{
	if (NOT Use_Greedy_Heuristic) {
		return smith_lee_liebman (pts, bsdp);
	}
	else {
		return greedy_heuristic	 (pts, bsdp);
	}
}

/*
 * Test if a new value of LP leads to a pruning of the eq-point.
 * If not then update LP.
 * First move the new point slightly to the left in order to
 * avoid (too many) numerical difficulties.
 */

	static
	bool
test_and_save_LP (

struct einfo *	eip,	/* IN - global EFST info */
struct eqp_t *	eqpi,	/* IN - first eq-point */
struct eqp_t *	eqpj,	/* IN - second eq-point */
struct eqp_t *	eqpk,	/* IN/OUT - new eq-point */
struct point *	CLP	/* IN - new proposed LP point */
)
{
struct point NLP;
struct point PLP;

	/* First check that the new point is between eqpk -> LP and eqpi -> E */
	if ((CLP -> x EQ eqpi -> E.x) AND (CLP -> y EQ eqpi -> E.y)) return TRUE;
	if (right_turn(&(eqpk -> E), &(eqpi -> E),  CLP))	     return TRUE;
	if (left_turn (&(eqpk -> E), &(eqpk -> LP), CLP))	     return TRUE;

	memset (&PLP, 0, sizeof (PLP));
	memset (&NLP, 0, sizeof (NLP));

	/* Now move it a little closer to eqpj -> E, just to be on the safe side */
	PLP.x = CLP -> x + (eqpj -> E.x - CLP -> x) * eip -> eps * ((double) eqpk -> S);
	PLP.y = CLP -> y + (eqpj -> E.y - CLP -> y) * eip -> eps * ((double) eqpk -> S);
	project_point(&(eqpk -> E), &(eqpk -> DC), &PLP, &NLP);

	/* We pass eqpk -> LP? */
	if (left_turn (&(eqpk -> E), &(eqpk -> LP), &NLP))	     return TRUE;

	/* Now make the actual testing */
	if (right_turn(&(eqpk -> E), &(eqpk -> RP), &NLP))	     return FALSE;
	eqpk -> LP = NLP;
	return TRUE;
}

/*
 * Test if a new value of RP leads to a pruning of the eq-point.
 * If not then update RP.
 * First move the new point slightly to the right in order to
 * avoid (too many) numerical difficulties.
 */

	static
	bool
test_and_save_RP (

struct einfo *	eip,	/* IN - global EFST info */
struct eqp_t *	eqpi,	/* IN - first eq-point */
struct eqp_t *	eqpj,	/* IN - second eq-point */
struct eqp_t *	eqpk,	/* IN/OUT - new eq-point */
struct point *	CRP	/* IN - new proposed RP point */
)
{
struct point NRP;
struct point PRP;

	/* First check that the new point is between eqpj -> E and eqpk -> RP */
	if ((CRP -> x EQ eqpj -> E.x) AND (CRP -> y EQ eqpj -> E.y)) return TRUE;
	if (left_turn (&(eqpk -> E), &(eqpj -> E),  CRP))	     return TRUE;
	if (right_turn(&(eqpk -> E), &(eqpk -> RP), CRP))	     return TRUE;

	memset (&PRP, 0, sizeof (PRP));
	memset (&NRP, 0, sizeof (NRP));

	/* Now move it a little closer to eqpi -> E, just to be on the safe side */
	PRP.x = CRP -> x + (eqpi -> E.x - CRP -> x) * eip -> eps * ((double) eqpk -> S);
	PRP.y = CRP -> y + (eqpi -> E.y - CRP -> y) * eip -> eps * ((double) eqpk -> S);
	project_point(&(eqpk -> E), &(eqpk -> DC), &PRP, &NRP);

	/* We pass eqpk -> RP? */
	if (right_turn (&(eqpk -> E), &(eqpk -> RP), &NRP))	     return TRUE;

	/* Now make the actual testing */
	if (left_turn(&(eqpk -> E), &(eqpk -> LP), &NRP))	     return FALSE;
	eqpk -> RP = NRP;
	return TRUE;
}

/*
 * Flag all terminals involved in eq-point k
 */

	static
	void
set_member_arr (

struct einfo * eip,  /* IN/OUT - global EFST info */
struct eqp_t * eqp,  /* IN - eq-point */
bool flag	     /* IN - value of flag */
)
{
	eterm_t *p    = eqp -> Z;
	eterm_t *endp = p + eqp -> S;

	while (p < endp) {
		int t = *p++;
		eip -> MEMB [t] = flag;
	}
}

/*