doprnt.c

Algorithms for Image Processing and Computer Vision Source Code

			for (n = fieldsz; n < dprec; n++)
				PUTC('0');

			/* the string or number proper */
			n = size;
			if (fp->_cnt - n >= 0 && (fp->_flag & _IOLBF) == 0) {
				fp->_cnt -= n;
				bcopy(t, (char *)fp->_ptr, n);
				fp->_ptr += n;
			} else
				while (--n >= 0)
					PUTC(*t++);
			/* trailing f.p. zeroes */
			while (--fpprec >= 0)
				PUTC('0');
			/* left-adjusting padding (always blank) */
			if (flags & LADJUST)
				for (n = realsz; n < width; n++)
					PUTC(' ');
			/* finally, adjust cnt */
			cnt += width > realsz ? width : realsz;
			break;
		case '\0':	/* "%?" prints ?, unless ? is NULL */
			return (cnt);
		default:
			PUTC((char)*fmt);
			cnt++;
		}
	}
	/* NOTREACHED */
}

static long double pten[] =
{
    1e1L, 1e2L, 1e4L, 1e8L, 1e16L, 1e32L, 1e64L, 1e128L, 1e256L,
    1e512L, 1e1024L, 1e2048L, 1e4096L
};

static long double ptenneg[] =
{
    1e-1L, 1e-2L, 1e-4L, 1e-8L, 1e-16L, 1e-32L, 1e-64L, 1e-128L, 1e-256L,
    1e-512L, 1e-1024L, 1e-2048L, 1e-4096L
};

#define MAXP 4096
#define NP   12
#define P    (4294967296.0L * 4294967296.0L * 2.0L)   /* 2^65 */
static long double INVPREC = P;
static long double PREC = 1.0L/P;
#undef P
/*
 * Defining FAST_LDOUBLE_CONVERSION results in a little bit faster
 * version, which might be less accurate (about 1 bit) for long
 * double. For 'normal' double it doesn't matter.
 */
/* #define FAST_LDOUBLE_CONVERSION */

static int
cvtl(long double number, int prec, int flags, char *signp, u_char fmtch,
     char *startp, char *endp)
{
	register char *p, *t;
	long double fract;
	int dotrim, expcnt, gformat;
	long double integer, tmp;

#ifdef __GO32__
	if ((expcnt = isspeciall(number, startp)))
		return(expcnt);
#endif

	dotrim = expcnt = gformat = 0;
	/* fract = modfl(number, &integer); */
	integer = number;

	/* get an extra slot for rounding. */
	t = ++startp;

	p = endp - 1;
	if (integer) {
		int i, lp=NP, pt=MAXP;
#ifndef FAST_LDOUBLE_CONVERSION
		long double oint = integer, dd=1.0L;
#endif
		if (integer > INVPREC) {
			integer *= PREC;
			while(lp >= 0) {
				if (integer >= pten[lp]) {
					expcnt += pt;
					integer *= ptenneg[lp];
#ifndef FAST_LDOUBLE_CONVERSION
					dd *= pten[lp];
#endif
				}
				pt >>= 1;
				lp--;
			}
#ifndef FAST_LDOUBLE_CONVERSION
			integer = oint/dd;
#else
			integer *= INVPREC;
#endif
		}
		/*
		 * Do we really need this ?
		 */
		for (i = 0; i < expcnt; i++)
		    *p-- = '0';
	}
	number = integer;
	fract = modfl(number, &integer);
	/*
	 * get integer portion of number; put into the end of the buffer; the
	 * .01 is added for modf(356.0 / 10, &integer) returning .59999999...
	 */
	for (; integer; ++expcnt) {
		tmp = modfl(integer * 0.1L , &integer);
		*p-- = tochar((int)((tmp + .01L) * 10));
	}
	switch(fmtch) {
	case 'f':
		/* reverse integer into beginning of buffer */
		if (expcnt)
			for (; ++p < endp; *t++ = *p);
		else
			*t++ = '0';
		/*
		 * if precision required or alternate flag set, add in a
		 * decimal point.
		 */
		if (prec || flags&ALT)
			*t++ = decimal;
		/* if requires more precision and some fraction left */
		if (fract) {
			if (prec)
				do {
					fract = modfl(fract * 10.0L, &tmp);
					*t++ = tochar((int)tmp);
				} while (--prec && fract);
			if (fract)
				startp = roundl(fract, (int *)NULL, startp,
				    t - 1, (char)0, signp);
		}
		for (; prec--; *t++ = '0');
		break;
	case 'e':
	case 'E':
eformat:	if (expcnt) {
			*t++ = *++p;
			if (prec || flags&ALT)
				*t++ = decimal;
			/* if requires more precision and some integer left */
			for (; prec && ++p < endp; --prec)
				*t++ = *p;
			/*
			 * if done precision and more of the integer component,
			 * round using it; adjust fract so we don't re-round
			 * later.
			 */
			if (!prec && ++p < endp) {
				fract = 0;
				startp = roundl((long double)0.0L, &expcnt,
						startp, t - 1, *p, signp);
			}
			/* adjust expcnt for digit in front of decimal */
			--expcnt;
		}
		/* until first fractional digit, decrement exponent */
		else if (fract) {
			int lp=NP, pt=MAXP;
#ifndef FAST_LDOUBLE_CONVERSION
			long double ofract = fract, dd=1.0L;
#endif
			expcnt = -1;
			if (fract < PREC) {
				fract *= INVPREC;
				while(lp >= 0) {
					if (fract <= ptenneg[lp]) {
						expcnt -= pt;
						fract *= pten[lp];
#ifndef FAST_LDOUBLE_CONVERSION
						dd *= pten[lp];
#endif
					}
					pt >>= 1;
					lp--;
				}
#ifndef FAST_LDOUBLE_CONVERSION
				fract = ofract*dd;
#else
				fract *= PREC;
#endif
			}
			/* adjust expcnt for digit in front of decimal */
			for (/* expcnt = -1 */ ;; --expcnt) {
				fract = modfl(fract * 10.0L, &tmp);
				if (tmp)
					break;
			}
			*t++ = tochar((int)tmp);
			if (prec || flags&ALT)
				*t++ = decimal;
		}
		else {
			*t++ = '0';
			if (prec || flags&ALT)
				*t++ = decimal;
		}
		/* if requires more precision and some fraction left */
		if (fract) {
			if (prec)
				do {
					fract = modfl(fract * 10.0L, &tmp);
					*t++ = tochar((int)tmp);
				} while (--prec && fract);
			if (fract)
				startp = roundl(fract, &expcnt, startp,
				    t - 1, (char)0, signp);
		}
		/* if requires more precision */
		for (; prec--; *t++ = '0');

		/* unless alternate flag, trim any g/G format trailing 0's */
		if (gformat && !(flags&ALT)) {
			while (t > startp && *--t == '0');
			if (*t == decimal)
				--t;
			++t;
		}
		t = exponentl(t, expcnt, fmtch);
		break;
	case 'g':
	case 'G':
		/* a precision of 0 is treated as a precision of 1. */
		if (!prec)
			++prec;
		/*
		 * ``The style used depends on the value converted; style e
		 * will be used only if the exponent resulting from the
		 * conversion is less than -4 or greater than the precision.''
		 *	-- ANSI X3J11
		 */
		if (expcnt > prec || (!expcnt && fract && fract < .0001)) {
			/*
			 * g/G format counts "significant digits, not digits of
			 * precision; for the e/E format, this just causes an
			 * off-by-one problem, i.e. g/G considers the digit
			 * before the decimal point significant and e/E doesn't
			 * count it as precision.
			 */
			--prec;
			fmtch -= 2;		/* G->E, g->e */
			gformat = 1;
			goto eformat;
		}
		/*
		 * reverse integer into beginning of buffer,
		 * note, decrement precision
		 */
		if (expcnt)
			for (; ++p < endp; *t++ = *p, --prec);
		else
			*t++ = '0';
		/*
		 * if precision required or alternate flag set, add in a
		 * decimal point.  If no digits yet, add in leading 0.
		 */
		if (prec || flags&ALT) {
			dotrim = 1;
			*t++ = decimal;
		}
		else
			dotrim = 0;
		/* if requires more precision and some fraction left */
		while (prec && fract) {
			fract = modfl(fract * 10.0L, &tmp);
			*t++ = tochar((int)tmp);
			prec--;
		}
		if (fract)
			startp = roundl(fract, (int *)NULL, startp, t - 1,
					(char)0, signp);
		/* alternate format, adds 0's for precision, else trim 0's */
		if (flags&ALT)
			for (; prec--; *t++ = '0');
		else if (dotrim) {
			while (t > startp && *--t == '0');
			if (*t != decimal)
				++t;
		}
	}
	return(t - startp);
}

static char *
roundl(long double fract, int *exp, char *start, char *end, char ch,
       char *signp)
{
	long double tmp;

	if (fract)
	{
		if (fract == 0.5L)
		{
		    char *e = end;
		    if (*e == '.')
			e--;
		    if (*e == '0' || *e == '2' || *e == '4'
			|| *e == '6' || *e == '8')
		    {
			tmp = 3.0;
			goto start;
		    }
		}
		(void)modfl(fract * 10.0L, &tmp);
	}
	else
		tmp = todigit(ch);
start:
	if (tmp > 4)
		for (;; --end) {
			if (*end == decimal)
				--end;
			if (++*end <= '9')
				break;
			*end = '0';
			if (end == start) {
				if (exp) {	/* e/E; increment exponent */
					*end = '1';
					++*exp;
				}
				else {		/* f; add extra digit */
					*--end = '1';
					--start;
				}
				break;
			}
		}
	/* ``"%.3f", (double)-0.0004'' gives you a negative 0. */
	else if (*signp == '-')
		for (;; --end) {
			if (*end == decimal)
				--end;
			if (*end != '0')
				break;
			if (end == start)
				*signp = 0;
		}
	return(start);
}

static char *
exponentl(char *p, int exp, u_char fmtch)
{
	register char *t;
	char expbuf[MAXEXPLD];

	*p++ = fmtch;
	if (exp < 0) {
		exp = -exp;
		*p++ = '-';
	}
	else
		*p++ = '+';
	t = expbuf + MAXEXPLD;
	if (exp > 9) {
		do {
			*--t = tochar(exp % 10);
		} while ((exp /= 10) > 9);
		*--t = tochar(exp);
		for (; t < expbuf + MAXEXPLD; *p++ = *t++);
	}
	else {
		*p++ = '0';
		*p++ = tochar(exp);
	}
	return(p);
}

#ifdef hp300
isspecial(d, bufp, signp)
	double d;
	char *bufp, *signp;
{
	register struct IEEEdp {
		unsigned sign:1;
		unsigned exp:11;
		unsigned manh:20;
		unsigned manl:32;
	} *ip = (struct IEEEdp *)&d;

	if (ip->exp != 0x7ff)
		return(0);
	if (ip->manh || ip->manl)
		(void)strcpy(bufp, "NaN");
	else
		(void)strcpy(bufp, "Inf");
	return(3);
}
#endif

#ifdef __GO32__
static int
isspeciall(long double d, char *bufp)
{
	register struct IEEExp {
        	unsigned manl:32;
		unsigned manh:32;
		unsigned exp:15;
        	unsigned sign:1;
	} *ip = (struct IEEExp *)&d;

	if (ip->exp != 0x7fff)
		return(0);
	if ((ip->manh & 0x7fffffff) || ip->manl)
    	{
		strcpy(bufp, "NaN");
        	nan = 1; /* kludge: we don't need the sign,  it's not nice
		            but it should work */
    	}
	else
		(void)strcpy(bufp, "Inf");
    	return(3);
}
#endif