p1read.c

生成直角Steiner树的程序包

		}
	}

	ip			= NEWA (total, int);
	cip -> term_trees	= NEWA (nverts + 1, int *);
	ptrs			= NEWA (nverts, int *);

	for (i = 0; i < nverts; i++) {
		cip -> term_trees [i]	= ip;
		ptrs [i]		= ip;
		ip += counts [i];
	}
	cip -> term_trees [i] = ip;

	for (i = 0; i < nedges; i++) {
		vp1 = cip -> edge [i];
		vp2 = cip -> edge [i + 1];
		while (vp1 < vp2) {
			j = *vp1++;
			*(ptrs [j])++ = i;
		}
	}

	free ((char *) ptrs);
	free ((char *) counts);
}

/*
 * This routine merges the FST incompatibility info we read in with
 * the "basic incompatibility" info we generate internally from the
 * FSTs themselves.  The final result is saved as "cinfo.inc_edges".
 */

	static
	void
init_inc_edges (

struct cinfo *		cip,		/* IN - compatibility info */
int **			incompat,	/* IN - input incompatibility lists */
int *			icounts		/* IN - lengths of incompat lists */
)
{
int			i;
int			j;
int			k;
int			nedges;
int			total;
int *			ip1;
int *			ip2;
int *			ip3;
int *			ip4;
int *			ip5;
int *			flist;
int *			newcounts;
int **			inc_edges;
int **			basic_incompat;
bitmap_t *		edge_mask;

	nedges = cip -> num_edges;
	edge_mask = cip -> initial_edge_mask;

	inc_edges = NEWA (nedges + 1, int *);
	newcounts = NEWA (nedges, int);

	basic_incompat = compute_basic_incompat (cip);

	flist = NEWA (nedges, int);

	total = 0;
	for (i = 0; i < nedges; i++) {
		k = icounts [i];
		ip1 = flist;
		ip2 = incompat [i];
		sort_ints (ip2, k);
		ip3 = ip2 + k;
		ip4 = basic_incompat [i];
		ip5 = basic_incompat [i + 1];
		for (;;) {
			for (;;) {
				if (ip2 >= ip3) break;
				j = *ip2;
				if (BITON (edge_mask, j)) break;
				++ip2;
			}
			if (ip2 >= ip3) {
				while (ip4 < ip5) {
					*ip1++ = *ip4++;
				}
				break;
			}
			if (ip4 >= ip5) {
				while (ip2 < ip3) {
					j = *ip2++;
					if (NOT BITON (edge_mask, j)) continue;
					*ip1++ = j;
				}
				break;
			}
			k = *ip4;
			if (j <= k) {
				if (j EQ k) {
					++ip4;
				}
				*ip1++ = j;
				++ip2;
			}
			else {
				*ip1++ = k;
				++ip4;
			}
		}
		k = ip1 - flist;
		newcounts [i] = k;
		total += k;
		ip2 = NULL;
		if (k > 0) {
			ip2 = NEWA (k, int);
			for (j = 0; j < k; j++) {
				ip2 [j] = flist [j];
			}
		}
		inc_edges [i] = ip2;
	}
	free((char *) flist);

	free (basic_incompat [0]);
	free (basic_incompat);

	ip1 = NEWA (total, int);
	for (i = 0; i < nedges; i++) {
		ip2 = inc_edges [i];
		inc_edges [i] = ip1;

		if (ip2 EQ NULL) continue;

		k = newcounts [i];
		for (j = 0; j < k; j++) {
			*ip1++ = ip2 [j];
		}
		free ((char *) ip2);
	}
	inc_edges [i] = ip1;

	free ((char *) newcounts);

	if (ip1 - inc_edges [0] NE total) {
		fatal ("init_inc_edges: Bug 1.");
	}

	cip -> inc_edges = inc_edges;

	verify_symmetric (cip -> inc_edges, nedges);
}

/*
 * This routine computes for each FST, a list of those FSTs having
 * at least 2 vertices in common.
 */

	int **
compute_basic_incompat (

struct cinfo *		cip	/* IN/OUT - compatibility info. */
)
{
int			i;
int			j;
int			k;
int			t;
int			fs;
int			common;
int			nedges;
int			kmasks;
int			nmasks;
int			total;
bitmap_t *		fsmask;
bitmap_t *		edge_mask;
bitmap_t *		tmask;
int *			ep1;
int *			ep2;
int *			vp1;
int *			vp2;
int *			flist;
int **			incompat;
int *			counts;

	kmasks	  = cip -> num_vert_masks;
	nedges	  = cip -> num_edges;
	nmasks	  = cip -> num_edge_masks;
	edge_mask = cip -> initial_edge_mask;

	incompat = NEWA (nedges + 1, int *);
	counts	 = NEWA (nedges, int);
	flist	 = NEWA (nedges, int);

	for (i = 0; i < nedges; i++) {
		incompat [i] = NULL;
		counts [i]   = 0;
	}

	fsmask = NEWA (nmasks, bitmap_t);
	for (i = 0; i < nmasks; i++) {
		fsmask [i] = 0;
	}

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

	/* Compute the list of (lists of) basically incomatible FSTs... */
	total = 0;
	for (i = 0; i < nedges; i++) {
		incompat [i] = NULL;
		if (NOT BITON (edge_mask, i)) continue;
		/* Develop list of all FSTs adjacent to FST i... */
		k = 0;
		vp1 = cip -> edge [i];
		vp2 = cip -> edge [i + 1];
		while (vp1 < vp2) {
			t = *vp1++;
			SETBIT (tmask, t);
			ep1 = cip -> term_trees [t];
			ep2 = cip -> term_trees [t + 1];
			while (ep1 < ep2) {
				fs = *ep1++;
				if (BITON (fsmask, fs)) continue;
				if (NOT BITON (edge_mask, fs)) continue;
				if (fs EQ i) continue;
				SETBIT (fsmask, fs);
				flist [k++] = fs;
			}
		}
		ep1 = &flist [0];
		ep2 = &flist [k];
		k = 0;
		while (ep1 < ep2) {
			fs = *ep1++;
			CLRBIT (fsmask, fs);
			/* Count number of vertices in common. */
			common = 0;
			vp1 = cip -> edge [fs];
			vp2 = cip -> edge [fs + 1];
			while (vp1 < vp2) {
				j = *vp1++;
				if (BITON (tmask, j)) {
					++common;
				}
			}
			if (common >= 2) {
				/* Too many in common!  Retain... */
				flist [k++] = fs;
			}
		}
		counts [i] = k;
		total += k;
		if (k > 0) {
			/* Save off sorted list of incompatible FSTs. */
			sort_ints (flist, k);
			ep1 = NEWA (k, int);
			incompat [i] = ep1;
			ep2 = flist;
			for (j = 0; j < k; j++) {
				*ep1++ = *ep2++;
			}
		}
		vp1 = cip -> edge [i];
		vp2 = cip -> edge [i + 1];
		while (vp1 < vp2) {
			j = *vp1++;
			CLRBIT (tmask, j);
		}
	}
	free ((char *) tmask);
	free ((char *) fsmask);
	free ((char *) flist);

	/* Now allocate and copy into contiguous memory... */
	ep1 = NEWA (total, int);
	for (i = 0; i < nedges; i++) {
		ep2 = incompat [i];
		incompat [i] = ep1;
		if (ep2 EQ NULL) continue;
		k = counts [i];
		for (j = 0; j < k; j++) {
			*ep1++ = ep2 [j];
		}
		free ((char *) ep2);
	}
	incompat [i] = ep1;

	free ((char *) counts);

	if (ep1 - incompat [0] NE total) {
		fatal ("compute_basic_incompat: Bug 1.");
	}

	return (incompat);
}

/*
 * Verify that the given "incompatibility" relation is symmetric.
 */

	static
	void
verify_symmetric (

int **		incompat,		/* IN - incompatibility lists */
int		nedges			/* IN - number of FSTs */
)
{
int		i;
int		j;
int **		begp;
int **		endp;

	begp = NEWA (nedges, int *);
	endp = NEWA (nedges, int *);

	for (i = 0; i < nedges; i++) {
		begp [i] = incompat [i];
		endp [i] = incompat [i + 1];
	}

	for (i = 0; i < nedges; i++) {
		for (;;) {
			if (begp [i] >= endp [i]) break;
			j = begp [i] [0];
			if (j > i) break;
			++(begp [i]);
			if ((begp [j] >= endp [j]) OR
			    (begp [j] [0] NE i)) {
				/* Incompatibility array not symmetric! */
				fatal ("verify_symmetric: Bug 1.");
			}
			++(begp [j]);
		}
	}

	for (i = 0; i < nedges; i++) {
		if (begp [i] NE endp [i]) {
			/* Incompatibility array not symmetric! */
			fatal ("verify_symmetric: Bug 2.");
		}
	}

	free ((char *) endp);
	free ((char *) begp);
}

/*
 * This routine determines the version number of the input data.
 */

	static
	int
get_version (void)

{
int		c;
int		n;
int		version;

	/* Strip any white space... */
	do {
		c = getchar ();
	} while ((c EQ ' ') OR
		 (c EQ '\t') OR
		 (c EQ '\n') OR
		 (c EQ '\f'));

	if (c < 0) {
		fprintf (stderr, "Unexpected EOF!\n");
		exit (1);
	}
	if (('0' <= c) AND (c <= '9')) {
		/* No version number input -- version 0. */
		ungetc (c, stdin);
		return (P1IO_VERSION_0);
	}

	if (c NE 'V') {
		fprintf (stderr, "Bad data version number!\n");
		exit (1);
	}
	version = -1;
	n = scanf ("%d", &version);
	if (n NE 1) {
		fprintf (stderr, "Data version number not found!\n");
		exit (1);
	}
	if (version <= P1IO_VERSION_0) {
		fprintf (stderr, "Corrupt version number!\n");
		exit (1);
	}
	if (version > CURRENT_P1IO_VERSION) {
		fprintf (stderr, "Version number out of range!\n");
		exit (1);
	}

	return (version);
}

/*
 * This routine reads a single decimal number from stdin.
 */

	static
	int
get_d (void)

{
int		i, n;

	n = scanf (" %d", &i);
	if (n NE 1) {
		fprintf (stderr, "Unexpected EOF!\n");
		exit (1);
	}

	return (i);
}


/*
 * This routine reads in a 32-bit bitmap element as a hex number
 * from stdin.
 */

	static
	bitmap_t
get_x (void)

{
bitmap_t	i;
int		n;

	n = scanf (" %lx", &i);
	if (n NE 1) {
		fprintf (stderr, "Unexpected EOF!\n");
		exit (1);
	}

	return (i);
}

/*
 * This routine reads in an entire line as a string.
 */

	static
	void
get_line (

char *		buf,		/* OUT - text of line */
int		buflen		/* IN - length of buffer */
)
{
int		n;
char *		p;

	p = fgets (buf, buflen, stdin);
	if (p EQ NULL) {
		fprintf (stderr, "Unexpected EOF!\n");
		exit (1);
	}
	n = strlen (buf);
	if (n > 0) {
		buf [n - 1] = '\0';	/* Discard newline. */
	}
}

/*
 * This routine reads a decimal floating point number.
 */

	static
	double
get_dec_double (void)

{
double		x;

	if (scanf (" %lg", &x) NE 1) {
		fprintf (stderr, "Expected floating point number.\n");
		exit (1);
	}

	return (x);
}

/*
 * This routine reads a HEXIDECIMAL floating point number.
 */

	static
	double
get_hex_double (void)

{
char		buf1 [128];

	if (scanf (" %s", buf1) NE 1) {
		fprintf (stderr, "Unexpected EOF!\n");
		exit (1);
	}

	return (hex_to_double (buf1));
}

/*
 * This routine skips the remainder of the current line.
 */

	static
	void
skip (void)

{
int		c;

	for (;;) {
		c = getchar ();
		if (c < 0) break;
		if (c EQ '\n') break;
	}
}

/*
 * This routine converts the printable ASCII hex format back into
 * a floating point number.
 */

	static
	double
hex_to_double (

char *		s		/* IN - the hex ASCII string to decode */
)
{
char		c;
char *		p1;
char *		p2;
int		digit;
int		expon;
int		esign;
double		msign;
double		mant;
double		value;

	/* Strip leading white space */
	for (;;) {
		c = *s++;
		if (c EQ ' ') continue;
		if (c EQ '\n') continue;
		if (c EQ '\t') continue;
		if (c EQ '\f') continue;
		if (c NE '\v') break;
	}

	msign = 1.0;
	if (c EQ '-') {
		msign = -1.0;
		c = *s++;
	}

	mant = 0.0;

	if (c NE '.') return (msign * mant);

	/* find end of mantissa */
	p1 = s;
	for (;;) {
		c = *s++;
		digit = hexdig (c);
		if (digit < 0) break;
	}

	/* rescan mantissa in reverse */
	for (p2 = s - 1; p2 > p1; ) {
		digit = hexdig (*--p2);
#if 0
		tracef ("	@	%d\n", digit);
#endif
		mant += ((double) digit);
		mant *= 0.0625;		/* divide by 16... */
	}

	expon = 0;
	if ((c EQ 'x') OR (c EQ 'X')) {
		c = *s++;
		esign = 1;
		if (c EQ '-') {
			esign = -1;
			c = *s++;
		}
		for (;;) {
			digit = hexdig (c);
			if (digit < 0) break;
			expon = (expon << 4) + digit;
			c = *s++;
		}
		expon *= esign;
	}

	value = msign * ldexp (mant, expon);

	return (value);
}

/*
 * Routine to identify and convert hexidecimal ASCII digits.
 */

	static
	int
hexdig (

char		c
)
{
	switch (c) {

	case '0':		return (0);
	case '1':		return (1);
	case '2':		return (2);
	case '3':		return (3);
	case '4':		return (4);
	case '5':		return (5);
	case '6':		return (6);
	case '7':		return (7);
	case '8':		return (8);
	case '9':		return (9);
	case 'A': case 'a':	return (10);
	case 'B': case 'b':	return (11);
	case 'C': case 'c':	return (12);
	case 'D': case 'd':	return (13);
	case 'E': case 'e':	return (14);
	case 'F': case 'f':	return (15);

	}

	return (-1);
}

/*
 * This routine turns off the "tmap" bit for all but the first
 * terminal in each duplicate terminal group.  This effectively
 * removes the duplicated terminals from the problem.
 */

	static
	void
remove_duplicates (

int		ndg,		/* IN - number of duplicate terminal groups */
int **		list,		/* IN - list of duplicate terminal groups */
bitmap_t *	tmap		/* IN/OUT - valid subset of terminals */
)
{
int		i;
int *		ip1;
int *		ip2;

	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) {
			CLRBIT (tmap, *ip1);
		}
	}
}

/*
 * This routine frees up the phase 1 data that was read in by
 * read_phase_1_data.
 */

	void
free_phase_1_data (

struct cinfo *		cip		/* IN - all phase 1 data */
)
{
int			i;
int			nedges;
struct full_set *	fsp;

	nedges	= cip -> num_edges;

	free ((char *) (cip -> edge [0]));
	free ((char *) (cip -> edge));
	free ((char *) (cip -> edge_size));
	free ((char *) (cip -> cost));
	free ((char *) (cip -> tflag));

	if (cip -> full_trees NE NULL) {
		for (i = 0; i < nedges; i++) {
			fsp = cip -> full_trees [i];
			free ((char *) (fsp -> terminals));
			free ((char *) (fsp -> steiners));
			free ((char *) (fsp -> edges));
			free ((char *) fsp);
		}
		free ((char *) (cip -> full_trees));
	}

	if (cip -> pts NE NULL) {
		free ((char *) (cip -> pts));
	}
	free ((char *) (cip -> initial_vert_mask));
	free ((char *) (cip -> initial_edge_mask));
	free ((char *) (cip -> required_edges));

	if (cip -> term_trees NE NULL) {
		if (cip -> term_trees [0] NE NULL) {
			free ((char *) (cip -> term_trees [0]));
		}
		free ((char *) (cip -> term_trees));
	}

	if (cip -> inc_edges NE NULL) {
		if (cip -> inc_edges [0] NE NULL) {
			free ((char *) (cip -> inc_edges [0]));
		}
		free ((char *) (cip -> inc_edges));
	}

	if (cip -> description NE NULL) {
		free ((char *) (cip -> description));
	}
}