io.c

生成直角Steiner树的程序包

double		mantissa,	/* IN - scaled mantissa value */
int		expon,		/* IN - exponent (power of 10) */
int		nsig		/* IN - num significant mantissa digits */
)
{
struct numlist *	p;

	p = NEW (struct numlist);

	p -> mantissa	= mantissa;
	p -> expon	= expon;
	p -> nsig	= (nsig > 0) ? nsig : 1;
	p -> next	= NULL;

	return (p);
}

/*
 * Compute a common scale value for all of the numbers in the given list.
 * In most cases, we can obtain an "optimum" scale value -- wherein every
 * number in the list will have an integer value that can be represented
 * exactly.  In cases where there are many significant digits, or a wide
 * span of exponent values, we must compromise.
 */

	int
compute_scaling_factor (

struct numlist *	p	/* IN - list of numbers to scale */
)
{
int		k;
int		min_expon;
int		max_expon;

	for (;;) {
		if (p EQ NULL) {
			/* No non-zero numbers -- don't scale at all! */
			return (0);
		}
		if (p -> mantissa NE 0.0) break;
		p = p -> next;
	}

	min_expon = p -> expon;
	max_expon = p -> expon + p -> nsig;

	for (p = p -> next; p NE NULL; p = p -> next) {
		if (p -> mantissa EQ 0) continue;
		if (p -> expon < min_expon) {
			min_expon = p -> expon;
		}
		k = p -> expon + p -> nsig;
		if (k > max_expon) {
			max_expon = k;
		}
	}

	/* Let D be a significant digit from the list, E be its		*/
	/* exponent so that the digit has value D * 10**E.  We now know	*/
	/* the smallest and largest of the E values.  If the spread is	*/
	/* MAX_INT_DIGITS or fewer, we win!  (Use the smallest E.)	*/
	/* Otherwise, try to split the difference.			*/

	max_expon -= MAX_INT_DIGITS;

	if (max_expon > min_expon) {
		/* Must compromise. */
		min_expon = (max_expon + min_expon + 1) / 2;
	}

	/* Our internal representation is such that dividing it by	*/
	/* 10**-min_expon yields external values.  (Note that min_expon	*/
	/* is typically negative here, in which case we would really be	*/
	/* multiplying by 10**(-min_expon).)				*/

	return (-min_expon);
}

/*
 * Set the scaling info properly for the given scale factor.
 */

	void
set_scale_info (

struct scale_info *	sip,		/* OUT - scale_info to initialize */
int			scale_factor	/* IN - scale factor to establish */
)
{
int		i;
double		pow10;
double		p;

	memset (sip, 0, sizeof (*sip));

	sip -> scale	= scale_factor;

	i = scale_factor;
	if (i < 0) {
		i = - i;
	}

	/* Compute p = 10 ** i. */
	pow10 = 1.0;
	p = 10.0;
	while (i > 0) {
		if ((i & 1) NE 0) {
			pow10 *= p;
		}
		p *= p;
		i >>= 1;
	}

	if (scale_factor < 0) {
		sip -> scale_mul = pow10;
		sip -> scale_div = 1.0;
	}
	else {
		sip -> scale_mul = 1.0;
		sip -> scale_div = pow10;
	}
}

/*
 * This routine converts the coordinate in the given numlist structure
 * to an internal coordinate value, scaled as specified by the scaling_factor
 * parameter.  The given numlist structure is freed.
 */

	coord_t
read_numlist (

struct numlist *	nlp,		/* IN - string list structure */
struct scale_info *	sip		/* IN - problem scaling info */
)
{
int		i;
int		target_expon;
int		expon;
coord_t		value;
double		pow10;
double		p;

	target_expon = - sip -> scale;

	value = nlp -> mantissa;
	expon = nlp -> expon;

	i = expon - target_expon;
	if (i < 0) {
		i = - i;
	}
	pow10 = 1.0;
	p = 10.0;
	while (i > 0) {
		if ((i & 1) NE 0) {
			pow10 *= p;
		}
		p *= p;
		i >>= 1;
	}
	if (expon > target_expon) {
		value *= pow10;
	}
	else if (expon < target_expon) {
		value /= pow10;
	}

	free ((char *) nlp);

	return (value);
}

/*
 * This routine gets the next character from stdin.  For the case of
 * interactive input, it makes sure that an EOF condition is
 * "permanent" -- that is, getc will not be called again after it
 * indicates EOF.
 */

	static
	int
next_char (

FILE *		fp		/* IN - input stream to read from */
)
{
int		c;

	if (feof (fp)) {
		/* Don't do "getc" again -- EOF is permanent! */
		return (EOF);
	}

	c = getc (fp);

	return (c);
}

/*
 * This routine converts the given internal-form coordinate to a
 * printable ASCII string in external form.  The internal form is an
 * integer (to eliminate numeric problems), but the external data
 * may actually involve decimal fractions.  This routine therefore
 * properly scales the coordinate during conversion to printable form.
 */

	void
coord_to_string (

char *			buf,	/* OUT - buffer to put ASCII text into */
coord_t			coord,	/* IN - coordinate to convert */
struct scale_info *	sip	/* IN - problem scaling info */
)
{
	/* Just pretend its a distance -- distances certainly	*/
	/* do not have LESS precision than coordinates...	*/
	dist_to_string (buf, coord, sip);
}

/*
 * This routine converts the given internal-form distance to a
 * printable ASCII string in external form.  The internal form is an
 * integer (to eliminate numeric problems), but the external data
 * may actually involve decimal fractions.  This routine therefore
 * properly scales the distance during conversion to printable form.
 */

	void
dist_to_string (

char *			buf,	/* OUT - buffer to put ASCII text into */
dist_t			dist,	/* IN - distance to convert */
struct scale_info *	sip	/* IN - problem scaling info */
)
{
int		i;
int		digit;
int		mindig;
double		ipart;
double		fpart;
double		div10;
char *		p;
char *		endp;
char		tmp [256];

	if (dist < 0.0) {
		dist = -dist;
		*buf++ = '-';
	}
	if (dist EQ 0) {
		*buf++ = '0';
		*buf++ = '\0';
		return;
	}

	mindig = sip -> min_precision;

	memset (tmp, '0', sizeof (tmp));
	p = &tmp [128];
	endp = p;

	ipart = floor (dist);
	while (ipart > 0) {
		div10 = floor (ipart / 10.0);
		digit = ipart - floor (10.0 * div10 + 0.5);
		*--p = digit + '0';
		ipart = div10;
	}
	if ((sip -> scale <= 0) AND (mindig EQ 0)) {
		/* The coordinates and edge costs are all integers	*/
		/* (both internally and externally so scaling is a	*/
		/* "no-op").  The metric is not Euclidean, so nothing	*/
		/* irrational appears.  Always output integers, and do	*/
		/* NOT insert a decimal point...			*/

		endp [-(sip -> scale)] = '\0';
		strcpy (buf, p);
		return;
	}

	if (sip -> scale >= 0) {
		if (p + sip -> scale > endp) {
			p = endp - sip -> scale;
		}

		for (i = 0; i <= sip -> scale; i++) {
			endp [1 - i] = endp [0 - i];
		}
		endp [1 - i] = '.';
		++endp;

		fpart = dist - floor (dist);
		while ((endp - p) < (mindig + 1)) {
			fpart *= 10.0;
			digit = (int) floor (fpart);
			*endp++ = digit + '0';
			fpart -= ((double) digit);
		}
	}
	else if (mindig EQ 0) {
		endp -= sip -> scale;
	}
	else {
		fpart = dist - floor (dist);
		for (i = sip -> scale; i < 0; i++) {
			fpart *= 10.0;
			digit = (int) floor (fpart);
			*endp++ = digit + '0';
			fpart -= ((double) digit);
			--mindig;
		}
		if (mindig > 0) {
			*endp++ = '.';
			while (mindig > 0) {
				fpart *= 10.0;
				digit = (int) floor (fpart);
				*endp++ = digit + '0';
				fpart -= ((double) digit);
				--mindig;
			}
		}
	}

	*endp = '\0';

	strcpy (buf, p);
}

/*
 * This routine sets up the various parameters needed for outputting
 * scaled coordinates.  We want to print coordinates/distances with
 * the minimum fixed precision whenever this gives the exact result.
 * Otherwise we will print the coordinates/distances with full
 * precision.
 *
 * The minimum fixed precision is exact ONLY when ALL of the following
 * are true:
 *
 *	- The metric must not be EUCLIDEAN, since distances become
 *	  irrational even if the coordinates are finite precision.
 *	  Coordinates of Euclidean Steiner points are also irrational.
 *
 *	- The SCALED vertex coordinates must all be integral.
 *
 *	- The SCALED hyperedge costs must all be integral.
 *
 * Note: we actually check the scaled data for integrality because there
 * are some old FST generators that do not implement scaling!  Such
 * data always contain a scale factor of zero and non-integral coordinates
 * and distances.
 */

	void
init_output_conversion (

struct pset *		pts,		/* IN - list of terminals */
int			metric,		/* IN - the metric */
struct scale_info *	sip		/* IN/OUT - problem scaling info */
)
{
int			i;
struct point *		p;
bool			integral;
double			c;

	integral = TRUE;

	if (pts NE NULL) {
		p = &(pts -> a [0]);
		for (i = 0; i < pts -> n; i++, p++) {
			c = floor ((double) (p -> x));
			if (c NE p -> x) {
				integral = FALSE;
				break;
			}
			c = floor ((double) (p -> y));
			if (c NE p -> y) {
				integral = FALSE;
				break;
			}
		}
	}

	if (integral AND (metric NE EUCLIDEAN)) {
		sip -> min_precision = 0;
	}
	else {
		sip -> min_precision = DBL_DIG + 1;
	}
}